Automation of absolute protein-ligand binding free energy calculations for docking refinement and compound evaluation

Assessment of the Topology and Oligomerisation States of Coiled Coils Using Metadynamics with Conformational Restraints

Coiled-coil proteins provide an excellent scaffold for multistate de novo protein design due to their established sequence-to-structure relationships and ability to switch conformations in response to external stimuli, such as changes in pH or temperature. However, the computational design of multistate coiled-coil protein assemblies is challenging, as it requires accurate estimates of the free energy differences between multiple alternative coiled-coil conformations. Here, we demonstrate how this challenge can be tackled using metadynamics simulations with orientational, positional and conformational restraints. We show that, even for subtle sequence variations, our protocol can predict the preferred topology of coiled-coil dimers and trimers, the preferred oligomerization states of coiled-coil dimers, trimers, and tetramers, as well as the switching behavior of a pH-dependent multistate system. Our approach provides a method for predicting the stability of coiled-coil designs and offers a new framework for computing binding free energies in protein-protein and multiprotein complexes.

FEP-SPell-ABFE: An Open-Source Automated Alchemical Absolute Binding Free-Energy Calculation Workflow for Drug Discovery

The binding affinity between a drug molecule and its target, measured by the absolute binding free energy (ABFE), is a crucial factor in the lead discovery phase of drug development. Recent research has highlighted the potential of in silico ABFE predictions to directly aid drug development by allowing for the ranking and prioritization of promising candidates. This work introduces an open-source Python workflow called FEP-SPell-ABFE, designed to automate ABFE calculations with minimal user involvement. The workflow requires only three key inputs: a receptor protein structure in PDB format, candidate ligands in SDF format, and a configuration file (config.yaml) that governs both the workflow and molecular dynamics simulation parameters. It produces a ranked list of ligands along with their binding free energies in the comma-separated values (CSV) format. The workflow leverages SLURM (Simple Linux Utility for Resource Management) for automating task execution and resource allocation across the modules. A usage example and several benchmark systems for validation are provided. The FEP-SPell-ABFE workflow, along with a practical example, is publicly accessible on GitHub at https://github.com/freeenergylab/FEP-SPell-ABFE, distributed under the MIT License.

Identification of Novel Fourth-Generation Allosteric Inhibitors Targeting Inactive State of EGFR T790M/L858R/C797S and T790M/L858R Mutations: A Combined Machine Learning and Molecular Dynamics Approach

Targeted therapy with an allosteric inhibitor (AIs) is an important area of research in patients with epidermal growth factor receptor (EGFR) mutations. Current treatment of nonsmall cell lung cancer patients with EGFR mutations using orthosteric inhibitors faces challenges like resistance and stopping over phosphorylation. Notably AIs have been introduced to overcome this resistance and increase inhibitory potency that binds to pockets other than the ATP-binding site (orthosteric site). Recently, fourth-generation AIs, EAI045, have been discovered to potently and selectively inhibit various EGFR mutations but limited antiproliferative effects in the absence of the antibody cetuximab. The purpose of this work is to identify nontoxic, potent small AIs through various screening pipelines and explore their molecular mechanism. In the discovery of AIs, structural similarity search, high-throughput virtual screening, and machine learning-guided QSAR modeling, several candidates were identified. Machine learning was employed to guide the QSAR model based on 2D descriptors and DFT-derived quantum chemical descriptors followed by a PCA reduction technique, which enabled the prediction of the biological activity (IC50) of screened drugs against various EGFR mutations such as T790M/L858R/C797S and T790M/L858R. In addition, multinanosecond (ns) and microsecond (μs) classical molecular dynamics (MD) simulations run on protein-ligand binding complex to check the stability of binding dynamics for T790M/L858R/C797S and T790M/L858R mutations with lower IC50 and higher docking score compounds. The molecular mechanics generalized Boltzmann surface area (MM/GBSA) calculation revealed that the five hit allosteric molecules for T790M/C797S/L858R and two for T790M/L858R mutations had a high binding affinity. The results were corroborated further by MM/GBSA employing the normal-mode analysis entropy method to perform additional screening. Furthermore, the compounds' efficacy was confirmed using path-dependent ligand unbinding free energy techniques such as Jarzynski averaged free energy profiles obtained from adaptive steered MD, relative residence time, and umbrella sampling simulations, which were compared to a reference inhibitor. However, path-independent alchemical approaches like streamlined alchemical free energy perturbation and binding free energy estimator 2 (BFEE2) were employed to validate the results and identify potent compounds. These findings pave the way to identification of novel potential fourth-generation AIs, which require further experimental validation.

SPONGE-FEP: An Automated Relative Binding Free Energy Calculation Accelerated by Selective Integrated Tempering Sampling

Computer-aided drug discovery (CADD) utilizes computational methods to accelerate the identification and optimization of potential drug candidates. Free energy perturbation (FEP) and thermodynamic integration (TI) play a critical role in predicting differences in protein binding affinities between drug molecules. Here, we implement SPONGE-FEP, which incorporates selective integrated tempering sampling (SITS) to enhance sampling efficiency and contains an automated workflow for relative binding free energy (RBFE) calculations. We first provide an overview of the workflow, which encompasses the generation of a perturbation map, alchemical free energy calculations, and cycle closure analysis. Two case studies were then performed to demonstrate the enhanced sampling of conformational states of ligands and proteins during the alchemical transformation process. The results show that the refined SITS method in SPONGE-FEP can significantly improve the sampling efficiency of rare events and the performance of RBFE predictions. Three series of comparative RBFE tests were conducted to demonstrate the accuracy of SPONGE-FEP, which is comparable to FEP+, using an average computation time of 4 h for a pair of ligands on an A100 GPU device.

Flexible framework of computing binding free energy using the energy representation theory of solution

Host-guest binding plays a crucial role in the functionality of various systems, and its efficiency is often quantified using the binding free energy, which represents the free-energy difference between the bound and dissociated states. Here, we propose a methodology to compute the binding free energy based on the energy representation (ER) theory of solution, which enables us to evaluate the free-energy difference between the systems of interest with the molecular dynamics (MD) simulations. Unlike the other free-energy methods, such as the Bennett acceptance ratio (BAR), the ER theory does not require the MD simulations for hypothetical intermediate states connecting the systems of interest, leading to reduced computational costs. By constructing the thermodynamic cycle of the binding process that is suitable for the ER theory, a robust calculation of the binding free energy is realized. We apply the present method to the self-association of N-methylacetamide in different solvents and the binding of aspirin to β-cyclodextrin (CD) in water. In the former case, the present method estimates that the binding free energy decreases as the solvent polarity decreases. This trend is consistent with the experimental finding. For the latter system, the binding free energies for the two representative CD-aspirin bound complexes, primary (P) and secondary (S) complexes, are estimated to be -5.2 ± 0.1 and -5.03 ± 0.09 kcal mol-1, respectively. These values are satisfactorily close to those from the BAR method [-4.2 ± 0.2 and -4.1 ± 0.2 kcal mol-1 for P and S, respectively]. Furthermore, the interaction-energy component analysis reveals that the van der Waals interaction between aspirin and CD dominantly contributes to the stabilization of the bound complexes, which is in harmony with the well-known binding mechanism in the CD systems.

ModBind, a Rapid Simulation-Based Predictor of Ligand Binding and Off-Rates

In rational drug discovery, both free energy of binding and the binding half-life (koff) are important factors in determining the efficacy of drugs. Numerous computational methods have been developed to predict these important properties, many of which rely on molecular dynamics (MD) simulations. While binding free-energy methods (thermodynamic equilibrium predictions) have been well validated and have demonstrated the ability to drive daily synthesis decisions in a commercial drug discovery setting, the prediction of koff (kinetics predictions) has had limited validation, and predictive methods have largely not been deployed in drug discovery settings. We developed ModBind, a novel method for MD simulation-based koff predictions. ModBind demonstrated similar accuracy to current state-of-the-art free-energy prediction methods. Additionally, ModBind performs ∼100 times faster than most available MD simulation-based free-energy or koff methods, allowing for widespread use by the molecular modeling community. While most free-energy methods rely on relative free-energy changes and are primarily useful for optimization of a congeneric series, our method requires no structural similarity between ligands, making ModBind an absolute predictor of koff. ModBind is thus a tool that can be used in virtual screening of diverse ligands, making it distinct from relative free-energy methods. We also discuss conditions that enable approximate prediction of ligand efficacy using ModBind and the limitations of this approach.

Essential Considerations for Free Energy Calculations of RNA-Small Molecule Complexes: Lessons from the Theophylline-Binding RNA Aptamer

Alchemical free energy calculations are widely used to predict the binding affinity of small molecule ligands to protein targets; however, the application of these methods to RNA targets has not been deeply explored. We systematically investigated how modeling decisions affect the performance of absolute binding free energy calculations for a relatively simple RNA model system: theophylline-binding RNA aptamer with theophylline and five analogs. The goal of this investigation was 2-fold: (1) understanding the performance levels we can expect from absolute free energy calculations for a simple RNA complex and (2) learning about practical modeling considerations that impact the success of RNA-binding predictions, which may be different from the best practices established for protein targets. We learned that magnesium ion (Mg2+) placement is a critical decision that impacts affinity predictions. When information regarding Mg2+ positions is lacking, implementing RNA backbone restraints is an alternative way of stabilizing the RNA structure that recapitulates prediction accuracy. Since mistakes in Mg2+ placement can be detrimental, omitting magnesium ions entirely and using RNA backbone restraints are attractive as a risk-mitigating approach. We found that predictions are sensitive to modeling experimental buffer conditions correctly, including salt type and ionic strength. We explored the effects of sampling in the alchemical protocol, choice of the ligand force field (GAFF2/OpenFF Sage), and water model (TIP3P/OPC) on predictions, which allowed us to give practical advice for the application of free energy methods to RNA targets. By capturing experimental buffer conditions and implementing RNA backbone restraints, we were able to compute binding affinities accurately (mean absolute error (MAE) = 2.2 kcal/mol, Pearson's correlation coefficient = 0.9, Kendall's τ = 0.7). We believe there is much to learn about how to apply free energy calculations for RNA targets and how to enhance their performance in prospective predictions. This study is an important first step for learning best practices and special considerations for RNA-ligand free energy calculations. Future studies will consider increasingly complicated ligands and diverse RNA systems and help the development of general protocols for therapeutically relevant RNA targets.

Comparative Assessment of Water Models in Protein-Glycan Interaction: Insights from Alchemical Free Energy Calculations and Molecular Dynamics Simulations

Accurate computational simulations of protein-glycan dynamics are crucial for a comprehensive understanding of critical biological mechanisms, including host-pathogen interactions, immune system defenses, and intercellular communication. The accuracy of these simulations, including molecular dynamics (MD) simulation and alchemical free energy calculations, critically relies on the appropriate parameters, including the water model, because of the extensive hydrogen bonding with glycan hydroxyl groups. However, a systematic evaluation of water models' accuracy in simulating protein-glycan interaction at the molecular level is still lacking. In this study, we used full atomistic MD simulations and alchemical absolute binding free energy (ABFE) calculations to investigate the performance of five distinct water models in six protein-glycan complex systems. We evaluated water models' impact on structural dynamics and binding affinity through over 5.8 μs of simulation time per system. Our results reveal that most protein-glycan complexes are stable in the overall structural dynamics regardless of the water model used, while some show obvious fluctuations with specific water models. More importantly, we discover that the stability of the binding motif's conformation is dependent on the water model chosen when its residues form weak hydrogen bonds with the glycan. The water model also influences the conformational stability of the glycan in its bound state according to density functional theory (DFT) calculations. Using alchemical ABFE calculations, we find that the OPC water model exhibits exceptional consistency with experimental binding affinity data, whereas commonly used models such as TIP3P are less accurate. The findings demonstrate how different water models affect protein-glycan interactions and the accuracy of binding affinity calculations, which is crucial in developing therapeutic strategies targeting these interactions.

Delving into Macrolide Binding Affinities and Associated Structural Modulations in Erythromycin Esterase C: Insights into the Venus Flytrap Mechanism

Since their inception in antibacterial therapy, macrolide-based antibiotics have significantly shaped the evolutionary pathways of pathogenic bacteria, driving them to develop diverse antimicrobial resistance (AMR) mechanisms. Among these, macrolide esterase, commonly referred to as erythromycin esterase, emerged as a critical defense mechanism, enabling bacteria to detoxify macrolides by hydrolyzing the macrolactone ring within the bacterial cell. In this study, we delve into the intricate interactions and conformational dynamics of erythromycin esterase C (EreC), a key member of the Ere enzyme family. We have focused on three FDA-approved and widely prescribed macrolides─erythromycin, clarithromycin, and azithromycin─by employing classical molecular dynamics, absolute binding free energy calculations, and 2D well-tempered metadynamics simulations to explore their interactions with EreC. To estimate the absolute binding free energies, we have used the recently developed and robust "Streamlined Alchemical Free Energy Perturbation (SAFEP)" protocol. The results from our molecular dynamics simulations and advanced analyses portrayed the crucial role of hydrophobic interactions within the macrolide binding cleft of EreC, along with the significant influence of the minor lobe in facilitating overall structural fluctuation. In silico alanine scanning identified top three hydrophobic residues, i.e., PHE248, MET333, and PHE344, responsible for macrolide binding inside that cleft. According to the free energy calculations, azithromycin and clarithromycin showed greater binding affinities toward EreC than the parent macrolide erythromycin. Moreover, 2D metadynamics simulations along with graph theory-based eigenvector centrality analyses revealed a metastable "semiopen" state during the hypothesized "active loop closure" of the EreC protein triggered by subtle conformational changes of an important histidine residue, HIS289, upon macrolide capture, drawing a fascinating parallel to the renowned "Venus flytrap" mechanism.

Konnektor: A Framework for Using Graph Theory to Plan Networks for Free Energy Calculations

Alchemical free energy campaigns can be planned using graph theory by building networks that contain nodes representing molecules that are connected by possible transformations as edges. We introduce Konnektor, an open-source Python package, for systematically planning, modifying, and analyzing free energy calculation networks. Konnektor is designed to aid in the drug discovery process by enabling users to easily setup free energy campaigns using complex graph manipulation methods. The package contains functions for network operations including concatenation of networks, deletion of transformations, and clustering of molecules along with a framework for combining these tools with existing network generation algorithms to enable the development of more complex methods for network generation. A comparison of the various network layout features offered is carried out using toy data sets. Additionally, Konnektor contains visualization and analysis tools, making the investigation of network features much simpler. Besides the content of the package, the paper also offers application examples, demonstrating how Konnektor can be used and how the different networks perform from a graph theory perspective. Konnektor is freely available via GitHub at https://github.com/OpenFreeEnergy/konnektor under the permissive MIT License.

In-silico assessment of novel peptidomimetics inhibitor targeting STAT3 and STAT4 N-terminal domain dimerization: A comprehensive study using molecular docking, molecular dynamics simulation, and binding free energy analysis

Dysregulation in Janus kinase-Signal Transducer and Activation of Transcription (JAK-STAT) pathway is closely linked to various cancer types. The N-terminal domain (NTD) of STAT proteins, upon dimerization, assumes a multifaceted role with remarkable adaptability in mediating interactions between proteins. Consequently, the strategic targeting of the N-terminal domain of STATs has emerged as a promising tactic for disrupting dimerization and impeding the translocation of STAT proteins. In this study, we have deployed an integrated in-silico methodology to rationally design Peptidomimetic foldamers as inhibitors of the N-terminal domains of STAT3 and STAT4, with the objective of disrupting protein dimerization. Consequently, we have judiciously designed a series of peptidomimetics that encompass β3-amino acids, bearing side chains that mimic the residues within interface II of the dimeric structures of the NTDs. Employing molecular docking techniques; we have assessed the binding affinity of these designed peptidomimetics toward both the NTDs. Furthermore, we have conducted an evaluation of the stability and conformational alterations within the docked complexes over an extensive Molecular Dynamics, subsequently computing the binding free energy utilizing MM/PBSA calculations. Our findings unequivocally demonstrate that the peptidomimetic foldamers we have devised (Peptide-A, Peptide-B, and Peptide-C) exhibit a propensity to bind to and impede the dimerization process of the NTDs of both STAT3 and STAT4. These outcomes serve to underscore the potential of these meticulously designed peptidomimetics as potential candidates meriting further exploration in the realm of cancer prevention and management.

Accurate protein-ligand binding free energy estimation using QM/MM on multi-conformers predicted from classical mining minima

Accurate prediction of binding free energy is crucial for the rational design of drug candidates and understanding protein-ligand interactions. To address this, we have developed four protocols that combine QM/MM calculations and the mining minima (M2) method, tested on 9 targets and 203 ligands. Our protocols carry out free energy processing with or without conformational search on the selected conformers obtained from M2 calculations, where their force field atomic charge parameters are substituted with those obtained from a QM/MM calculation. The method achieved a high Pearson's correlation coefficient (0.81) with experimental binding free energies across diverse targets, demonstrating its generality. Using a differential evolution algorithm with a universal scaling factor of 0.2, we achieved a low mean absolute error of 0.60 kcal mol-1. This performance surpasses many existing methods and is comparable to popular relative binding free energy techniques but at significantly lower computational cost.

Alchemical Transformations and Beyond: Recent Advances and Real-World Applications of Free Energy Calculations in Drug Discovery

Computational methods constitute efficient strategies for screening and optimizing potential drug molecules. A critical factor in this process is the binding affinity between candidate molecules and targets, quantified as binding free energy. Among various estimation methods, alchemical transformation methods stand out for their theoretical rigor. Despite challenges in force field accuracy and sampling efficiency, advancements in algorithms, software, and hardware have increased the application of free energy perturbation (FEP) calculations in the pharmaceutical industry. Here, we review the practical applications of FEP in drug discovery projects since 2018, covering both ligand-centric and residue-centric transformations. We show that relative binding free energy calculations have steadily achieved chemical accuracy in real-world applications. In addition, we discuss alternative physics-based simulation methods and the incorporation of deep learning into free energy calculations.

Optimal Molecular Design: Generative Active Learning Combining REINVENT with Precise Binding Free Energy Ranking Simulations

Active learning (AL) is a specific instance of sequential experimental design and uses machine learning to intelligently choose the next data point or batch of molecular structures to be evaluated. In this sense, it closely mimics the iterative design-make-test-analysis cycle of laboratory experiments to find optimized compounds for a given design task. Here, we describe an AL protocol which combines generative molecular AI, using REINVENT, and physics-based absolute binding free energy molecular dynamics simulation, using ESMACS, to discover new ligands for two different target proteins, 3CLpro and TNKS2. We have deployed our generative active learning (GAL) protocol on Frontier, the world's only exa-scale machine. We show that the protocol can find higher-scoring molecules compared to the baseline, a surrogate ML docking model for 3CLpro and compounds with experimentally determined binding affinities for TNKS2. The ligands found are also chemically diverse and occupy a different chemical space than the baseline. We vary the batch sizes that are put forward for free energy assessment in each GAL cycle to assess the impact on their efficiency on the GAL protocol and recommend their optimal values in different scenarios. Overall, we demonstrate a powerful capability of the combination of physics-based and AI methods which yields effective chemical space sampling at an unprecedented scale and is of immediate and direct relevance to modern, data-driven drug discovery.

Machine Learning Guided AQFEP: A Fast and Efficient Absolute Free Energy Perturbation Solution for Virtual Screening

Structure-based methods in drug discovery have become an integral part of the modern drug discovery process. The power of virtual screening lies in its ability to rapidly and cost-effectively explore enormous chemical spaces to select promising ligands for further experimental investigation. Relative free energy perturbation (RFEP) and similar methods are the gold standard for binding affinity prediction in drug discovery hit-to-lead and lead optimization phases, but have high computational cost and the requirement of a structural analog with a known activity. Without a reference molecule requirement, absolute FEP (AFEP) has, in theory, better accuracy for hit ID, but in practice, the slow throughput is not compatible with VS, where fast docking and unreliable scoring functions are still the standard. Here, we present an integrated workflow to virtually screen large and diverse chemical libraries efficiently, combining active learning with a physics-based scoring function based on a fast absolute free energy perturbation method. We validated the performance of the approach in the ranking of structurally related ligands, virtual screening hit rate enrichment, and active learning chemical space exploration; disclosing the largest reported collection of free energy simulations to date.

BAT2: an Open-Source Tool for Flexible, Automated, and Low Cost Absolute Binding Free Energy Calculations

Absolute binding free energy (ABFE) calculations with all-atom molecular dynamics (MD) have the potential to greatly reduce costs in the first stages of drug discovery. Here, we introduce BAT2, the new version of the Binding Affinity Tool (BAT.py), designed to combine full automation of ABFE calculations with high-performance MD simulations, making it a potential tool for virtual screening. We describe and test several changes and new features that were incorporated into the code, such as relative restraints between the protein and the ligand instead of using fixed dummy atoms, support for the OpenMM simulation engine, a merged approach to the application/release of restraints, support for cobinders and proteins with multiple chains, and many others. We also reduced the simulation times for each ABFE calculation, assessing the effect on the expected robustness and accuracy of the calculations.

Repurposing of antibacterial compounds for suppression of Mycobacterium tuberculosis dormancy reactivation by targeting resuscitation-promoting factors B

Tuberculosis infection has always been a global concern for public health, and the mortality rate has increased tremendously every year. The ability of the resuscitation Mycobacterium tuberculosis (Mtb) from the dormant state is one of the major reasons for the epidemic spread of tuberculosis infection, especially latent tuberculosis infection (LTBI). The element that encourages resuscitation, RpfB (resuscitation-promoting factors B), is mostly in charge of bringing Mtb out of slumber. This reason makes RpfB a promising target for developing tuberculosis drugs because of the effects of latent tuberculosis. Therefore, this work was executed using a computational three-level screening of the Selleckhem antibiotics database consisting of 462 antibiotics against the ligand binding region of the RpfB protein, followed by an estimation of binding free energy for ideal identification and confirmation of potential RpfB inhibitor. Subsequently, three antibiotic drug molecules, i.e., Amikacin hydrate (-66.87 kcal/mol), Isepamicin sulphate (-60.8 kcal/mol), and Bekanamycin (-46.89 kcal/mol), were selected on the basis of their binding free energy value for further computational studies in comparison to reference ligand, 4-benzoyl-2-nitrophenyl thiocyanate (NPT7). Based on the intermolecular interaction profiling, 200 ns molecular dynamic simulation (MD), post-simulation analysis and principal component analysis (PCA), the selected antibiotics showed substantial stability with the RpfB protein compared to the NPT7 inhibitor. Conclusively based on the computational results, the preferred drugs can be potent inhibitors of the RpfB protein, which can be further validated using in vivo research and in vitro enzyme inhibition to understand their therapeutic activity against tuberculosis infection.Communicated by Ramaswamy H. Sarma.

The molecular and network mechanisms of antilipidemic potential effects of Ganfule capsules in nonalcoholic fatty liver disease

Background

Non-alcoholic fatty liver disease (NAFLD) is a common liver disorder characterized by hepatic steatosis, inflammation and fibrosis. Ganfule (GFL), a traditional Chinese medicine, has demonstrated therapeutic potential in the treatment of NAFLD but the mechanisms involved are not fully understood.To evaluate the biochemical mechanisms of GFL in treating NAFLD by examining its effects on biological networks, key therapeutic targets, histopathological changes and clinical implications.

Methods

Chemical component screening, key target prediction, biological functional enrichment analysis, lipid profile localization analysis and complex network analysis were performed on GFL using multi-database mining, network analysis and molecular docking. An NAFLD rat model was then established and treated with different doses of GFL. Histopathological evaluation and western blotting were used to verify the expression levels of key target proteins in GFL-treated NAFLD rats.

Results

Network analysis analysis identified 12 core targets, 12 core active ingredients and 7 core Chinese medicinal herbs in GFL potentially involved in the treatment of NAFLD. Biological functional enrichment analysis revealed the involvement of lipid metabolism, apoptosis and intracellular signaling pathways. Molecular docking confirmed a strong affinity between GFL's core compounds and certain target proteins. Histopathological examination of an NAFLD rat model showed reduced hepatocellular steatosis after GFL treatment. Western blotting revealed significant downregulation of PPARA and PPARD protein expression and upregulation of PIK3CG and PRKACA protein expression in NAFLD rats treated with lower doses of GFL.

Conclusions

Our results suggest that GFL modulates key proteins involved in lipid metabolism and apoptosis pathways. GFL improved the histopathological features of NAFLD rats by regulating lipid metabolism as well as reducing hepatocyte apoptosis and hepatocellular steatosis. These findings offer insights into the biochemical mechanism of action of GFL and support its use in the treatment for NAFLD.

Automated Absolute Binding Free Energy Calculation Workflow for Drug Discovery

Absolute binding free energies play a crucial role in drug development, particularly as part of the lead discovery process. In recent work, we showed how in silico predictions directly could support drug development by ranking and recommending favorable ideas over unfavorable ones. Here, we demonstrate a Python workflow that enables the calculation of ABFEs with minimal manual input effort, such as the receptor PDB and ligand SDF files, and outputs a .tsv file containing the ranked ligands and their corresponding binding free energies. The implementation uses Snakemake to structure and control the execution of tasks, allowing for dynamic control of parameters and execution patterns. We provide an example of a benchmark system that demonstrates the effectiveness of the automated workflow.

Kartograf: A Geometrically Accurate Atom Mapper for Hybrid-Topology Relative Free Energy Calculations

Relative binding free energy (RBFE) calculations have emerged as a powerful tool that supports ligand optimization in drug discovery. Despite many successes, the use of RBFEs can often be limited by automation problems, in particular, the setup of such calculations. Atom mapping algorithms are an essential component in setting up automatic large-scale hybrid-topology RBFE calculation campaigns. Traditional algorithms typically employ a 2D subgraph isomorphism solver (SIS) in order to estimate the maximum common substructure. SIS-based approaches can be limited by time-intensive operations and issues with capturing geometry-linked chemical properties, potentially leading to suboptimal solutions. To overcome these limitations, we have developed Kartograf, a geometric-graph-based algorithm that uses primarily the 3D coordinates of atoms to find a mapping between two ligands. In free energy approaches, the ligand conformations are usually derived from docking or other previous modeling approaches, giving the coordinates a certain importance. By considering the spatial relationships between atoms related to the molecule coordinates, our algorithm bypasses the computationally complex subgraph matching of SIS-based approaches and reduces the problem to a much simpler bipartite graph matching problem. Moreover, Kartograf effectively circumvents typical mapping issues induced by molecule symmetry and stereoisomerism, making it a more robust approach for atom mapping from a geometric perspective. To validate our method, we calculated mappings with our novel approach using a diverse set of small molecules and used the mappings in relative hydration and binding free energy calculations. The comparison with two SIS-based algorithms showed that Kartograf offers a fast alternative approach. The code for Kartograf is freely available on GitHub (https://github.com/OpenFreeEnergy/kartograf). While developed for the OpenFE ecosystem, Kartograf can also be utilized as a standalone Python package.

Toward physics-based precision medicine: Exploiting protein dynamics to design new therapeutics and interpret variants

The goal of precision medicine is to utilize our knowledge of the molecular causes of disease to better diagnose and treat patients. However, there is a substantial mismatch between the small number of food and drug administration (FDA)-approved drugs and annotated coding variants compared to the needs of precision medicine. This review introduces the concept of physics-based precision medicine, a scalable framework that promises to improve our understanding of sequence-function relationships and accelerate drug discovery. We show that accounting for the ensemble of structures a protein adopts in solution with computer simulations overcomes many of the limitations imposed by assuming a single protein structure. We highlight studies of protein dynamics and recent methods for the analysis of structural ensembles. These studies demonstrate that differences in conformational distributions predict functional differences within protein families and between variants. Thanks to new computational tools that are providing unprecedented access to protein structural ensembles, this insight may enable accurate predictions of variant pathogenicity for entire libraries of variants. We further show that explicitly accounting for protein ensembles, with methods like alchemical free energy calculations or docking to Markov state models, can uncover novel lead compounds. To conclude, we demonstrate that cryptic pockets, or cavities absent in experimental structures, provide an avenue to target proteins that are currently considered undruggable. Taken together, our review provides a roadmap for the field of protein science to accelerate precision medicine.

Addressing Suboptimal Poses in Nonequilibrium Alchemical Calculations

Alchemical transformations can be used to quantitatively estimate absolute binding free energies at a reasonable computational cost. However, most of the approaches currently in use require knowledge of the correct (crystallographic) pose. In this paper, we present a combined Hamiltonian replica exchange nonequilibrium alchemical method that allows us to reliably calculate absolute binding free energies, even when starting from suboptimal initial binding poses. Performing a preliminary Hamiltonian replica exchange enhances the sampling of slow degrees of freedom of the ligand and the target, allowing the system to populate the correct binding pose when starting from an approximate docking pose. We apply the method on 6 ligands of the first bromodomain of the BRD4 bromodomain-containing protein. For each ligand, we start nonequilibrium alchemical transformations from both the crystallographic pose and the top-scoring docked pose that are often significantly different. We show that the method produces statistically equivalent binding free energies, making it a useful tool for computational drug discovery pipelines.

Tuning Potential Functions to Host-Guest Binding Data

Software to more rapidly and accurately predict protein-ligand binding affinities is of high interest for early-stage drug discovery, and physics-based methods are among the most widely used technologies for this purpose. The accuracy of these methods depends critically on the accuracy of the potential functions that they use. Potential functions are typically trained against a combination of quantum chemical and experimental data. However, although binding affinities are among the most important quantities to predict, experimental binding affinities have not to date been integrated into the experimental data set used to train potential functions. In recent years, the use of host-guest complexes as simple and tractable models of binding thermodynamics has gained popularity due to their small size and simplicity, relative to protein-ligand systems. Host-guest complexes can also avoid ambiguities that arise in protein-ligand systems such as uncertain protonation states. Thus, experimental host-guest binding data are an appealing additional data type to integrate into the experimental data set used to optimize potential functions. Here, we report the extension of the Open Force Field Evaluator framework to enable the systematic calculation of host-guest binding free energies and their gradients with respect to force field parameters, coupled with the curation of 126 host-guest complexes with available experimental binding free energies. As an initial application of this novel infrastructure, we optimized generalized Born (GB) cavity radii for the OBC2 GB implicit solvent model against experimental data for 36 host-guest systems. This refitting led to a dramatic improvement in accuracy for both the training set and a separate test set with 90 additional host-guest systems. The optimized radii also showed encouraging transferability from host-guest systems to 59 protein-ligand systems. However, the new radii are significantly smaller than the baseline radii and lead to excessively favorable hydration free energies (HFEs). Thus, users of the OBC2 GB model currently may choose between GB cavity radii that yield more accurate binding affinities and GB cavity radii that yield more accurate HFEs. We suspect that achieving good accuracy on both will require more far-reaching adjustments to the GB model. We note that binding free-energy calculations using the OBC2 model in OpenMM gain about a 10× speedup relative to corresponding explicit solvent calculations, suggesting a future role for implicit solvent absolute binding free-energy (ABFE) calculations in virtual compound screening. This study proves the principle of using host-guest systems to train potential functions that are transferrable to protein-ligand systems and provides an infrastructure that enables a range of applications.

Alchemical Free Energy Workflows for the Computation of Protein-Ligand Binding Affinities

Alchemical free energy methods can be used for the efficient computation of relative binding free energies during preclinical drug discovery stages. In recent years, this has been facilitated further by the implementation of workflows that enable non-experts to quickly and consistently set up the required simulations. Given the correct input structures, workflows handle the difficult aspects of setting up perturbations, including consistently defining the perturbable molecule, its atom mapping and topology generation, perturbation network generation, running of the simulations via different sampling methods, and analysis of the results. Different academic and commercial workflows are discussed, including FEW, FESetup, FEPrepare, CHARMM-GUI, Transformato, PMX, QLigFEP, TIES, ProFESSA, PyAutoFEP, BioSimSpace, FEP+, Flare, and Orion. These workflows differ in various aspects, such as mapping algorithms or enhanced sampling methods. Some workflows can accommodate more than one molecular dynamics (MD) engine and use external libraries for tasks. Differences between workflows can present advantages for different use cases, however a lack of interoperability of the workflows' components hinders systematic comparisons.

Shape-Based Virtual Screening of a Billion-Compound Library Identifies Mycobacterial Lipoamide Dehydrogenase Inhibitors

Lpd (lipoamide dehydrogenase) in Mycobacterium tuberculosis (Mtb) is required for virulence and is a genetically validated tuberculosis (TB) target. Numerous screens have been performed over the last decade, yet only two inhibitor series have been identified. Recent advances in large-scale virtual screening methods combined with make-on-demand compound libraries have shown the potential for finding novel hits. In this study, the Enamine REAL library consisting of ∼1.12 billion compounds was efficiently screened using the GPU Shape screen method against Mtb Lpd to find additional chemical matter that would expand on the known sulfonamide inhibitor series. We identified six new inhibitors with IC50 in the range of 5-100 μM. While these compounds remained chemically close to the already known sulfonamide series inhibitors, some diversity was found in the cores of the hits. The two most potent hits were further validated by one-step potency optimization to submicromolar levels. The co-crystal structure of optimized analogue TDI-13537 provided new insights into the potency determinants of the series.

Standard Binding Free-Energy Calculations: How Far Are We from Automation?

Recent success stories suggest that in silico protein-ligand binding free-energy calculations are approaching chemical accuracy. However, their widespread application remains limited by the extensive human intervention required, posing challenges for the neophyte. As such, it is critical to develop automated workflows for estimating protein-ligand binding affinities with minimum personal involvement. Key human efforts include setting up and tuning enhanced-sampling or alchemical-transformation algorithms as a preamble to computational binding free-energy estimations. Additionally, preparing input files, bookkeeping, and postprocessing represent nontrivial tasks. In this Perspective, we discuss recent progress in automating standard binding free-energy calculations, featuring the development of adaptive or parameter-free algorithms, standardization of binding free-energy calculation workflows, and the implementation of user-friendly software. We also assess the current state of automated standard binding free-energy calculations and evaluate the limitations of existing methods. Last, we outline the requirements for future algorithms and workflows to facilitate automated free-energy calculations for diverse protein-ligand complexes.

Landscape of In Silico Tools for Modeling Covalent Modification of Proteins: A Review on Computational Covalent Drug Discovery

Covalent drug discovery has been a challenging research area given the struggle of finding a sweet balance between selectivity and reactivity for these drugs, the lack of which often leads to off-target activities and hence undesirable side effects. However, there has been a resurgence in covalent drug design following the success of several covalent drugs such as boceprevir (2011), ibrutinib (2013), neratinib (2017), dacomitinib (2018), zanubrutinib (2019), and many others. Design of covalent drugs includes many crucial factors, where "evaluation of the binding affinity" and "a detailed mechanistic understanding on covalent inhibition" are at the top of the list. Well-defined experimental techniques are available to elucidate these factors; however, often they are expensive and/or time-consuming and hence not suitable for high throughput screens. Recent developments in in silico methods provide promise in this direction. In this report, we review a set of recent publications that focused on developing and/or implementing novel in silico techniques in "Computational Covalent Drug Discovery (CCDD)". We also discuss the advantages and disadvantages of these approaches along with what improvements are required to make it a great tool in medicinal chemistry in the near future.

Computation of Absolute Binding Free Energies for Noncovalent Inhibitors with SARS-CoV-2 Main Protease

Accurate, routine calculation of absolute binding free energies (ABFEs) for protein-ligand complexes remains a key goal of computer-aided drug design since it can enable screening and optimization of drug candidates. For development and testing of related methods, it is important to have high-quality datasets. To this end, from our own experimental studies, we have selected a set of 16 inhibitors of the SARS-CoV-2 main protease (Mpro) with structural diversity and well-distributed BFEs covering a 5 kcal/mol range. There is also minimal structural uncertainty since X-ray crystal structures have been deposited for 12 of the compounds. For methods testing, we report ABFE results from 2 μs molecular dynamics (MD) simulations using free energy perturbation (FEP) theory. The correlation of experimental and computed results is encouraging, with a Pearson's r2 of 0.58 and a Kendall τ of 0.24. The results indicate that current FEP-based ABFE calculations can be used for identification of active compounds (hits). While their accuracy for lead optimization is not yet sufficient, this activity remains addressable in separate lead series by relative BFE calculations.

Comparison of Receptor-Ligand Restraint Schemes for Alchemical Absolute Binding Free Energy Calculations

Alchemical absolute binding free energy calculations are of increasing interest in drug discovery. These calculations require restraints between the receptor and ligand to restrict their relative positions and, optionally, orientations. Boresch restraints are commonly used, but they must be carefully selected in order to sufficiently restrain the ligand and to avoid inherent instabilities. Applying multiple distance restraints between anchor points in the receptor and ligand provides an alternative framework without inherent instabilities which may provide convergence benefits by more strongly restricting the relative movements of the receptor and ligand. However, there is no simple method to calculate the free energy of releasing these restraints due to the coupling of the internal and external degrees of freedom of the receptor and ligand. Here, a method to rigorously calculate free energies of binding with multiple distance restraints by imposing intramolecular restraints on the anchor points is proposed. Absolute binding free energies for the human macrophage migration inhibitory factor/MIF180, system obtained using a variety of Boresch restraints and rigorous and nonrigorous implementations of multiple distance restraints are compared. It is shown that several multiple distance restraint schemes produce estimates in good agreement with Boresch restraints. In contrast, calculations without orientational restraints produce erroneously favorable free energies of binding by up to approximately 4 kcal mol-1. These approaches offer new options for the deployment of alchemical absolute binding free energy calculations.

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

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

Combining QM/MM Calculations with Classical Mining Minima to Predict Protein-Ligand Binding Free Energy

We developed an effective binding free energy prediction protocol which incorporates quantum mechanical/molecular mechanical (QM/MM) calculations to substitute the specified atomic charges of force fields with quantum-mechanically recalculated ones at a proposed pose using a mining minima approach with the VeraChem mining minima engine. We tested this protocol using seven well-known targets with 147 different ligands and compared it with classical mining minima and the most popular binding free energy (BFE) methods using different metrics. Our new protocol, dubbed Qcharge-VM2, yielded an overall Pearson correlation of 0.86, which was better than all the methods examined. Qcharge-VM2 performed significantly better than implicit solvent-based methods, such as MM-GBSA and MM-PBSA, but not as good as explicit water-based free energy perturbation methods, such as FEP+, in terms of root-mean-square error, RMSE (1.75 kcal/mol) and mean unsigned error, MUE (1.39 kcal/mol) on a limited set of targets. However, our protocol is substantially less computationally demanding compared with FEP+. The combined accuracy and efficiency of our method can be valuable in drug discovery campaigns.

Investigating the Solvent Effects on Binding Affinity of PAHs-ExBox4+ Complexes: An Alchemical Approach

Polycyclic aromatic hydrocarbons (PAHs) are polluting agents, produced naturally or artificially, widely dispersed in the environment and potentially carcinogenic and immunotoxic to humans and animals, mainly for marine life. Recently, a tetracationic box-shaped cyclophane (ExBox⁴⁺) was synthesized, fully characterized, and revealed to form host-guest complexes with PAHs in acetonitrile, demonstrating the potential ability for it to act as a PAHs scavenger. This work investigates, through Molecular Dynamics (MD) simulations, the binding affinity between different PAHs and ExBox⁴⁺ in different solvents: chloroform (nonpolar), acetonitrile (polar protic), and water (polar protic). An alchemical method of simultaneous decoupling-recoupling (SDR) was used and implemented in a newly developed Python program called GHOAT, which fully automates the calculation of binding free energies and invokes the AMBER 2020 simulation package. The results showed that the affinity between ExBox⁴⁺ and PAHs in water is much larger than in organic media, with free energies between -5 and -20 kcal/mol, being able to act as a PHAs scavenger with great potential for applications in environmental chemistry such as soil washing. The results also reveal a significant correlation with the experimental available ΔG values. The methodology employed presents itself as an important tool for the in silico determination of binding affinities, not only available for charged cyclophanes but also extensible to several other HG supramolecular systems in condensed media, aiding in the rational design of host-guest systems in a significant way.

Capturing the Recognition Dynamics of para-Sulfonato-calix[4]arenes by Cytochrome c: Toward a Quantitative Free Energy Assessment

Calix[n]arenes' selective recognition of protein surfaces covers a broad range of timely applications, from controlling protein assembly and crystallization to trapping partially disordered proteins. Here, the interaction of para-sulfonated calix-[4]-arenes with cytochrome c is investigated through all-atom, explicit water molecular dynamics simulations which allow characterization of two binding sites in quantitative agreement with experimental evidence. Free energy calculations based on the MM-PBSA and the attach-pull-release (APR) methods highlight key residues implicated in the recognition process and provide binding free energy results in quantitative agreement with isothermal titration calorimetry. Our study emphasizes the role of MD simulations to capture and describe the "walk" of sulfonated calix-[4]-arenes on the cytochrome c surface, with the arginine R13 as a pivotal interacting residue. Our MD investigation allows, through the quasi-harmonic multibasin (QHMB) method, probing an allosteric reinforcement of several per-residue interactions upon calixarene binding, which suggests a more complex mode of action of these supramolecular auxiliaries.

Evolutionary divergence in the conformational landscapes of tyrosine vs serine/threonine kinases

Inactive conformations of protein kinase catalytic domains where the DFG motif has a "DFG-out" orientation and the activation loop is folded present a druggable binding pocket that is targeted by FDA-approved 'type-II inhibitors' in the treatment of cancers. Tyrosine kinases (TKs) typically show strong binding affinity with a wide spectrum of type-II inhibitors while serine/threonine kinases (STKs) usually bind more weakly which we suggest here is due to differences in the folded to extended conformational equilibrium of the activation loop between TKs vs. STKs. To investigate this, we use sequence covariation analysis with a Potts Hamiltonian statistical energy model to guide absolute binding free-energy molecular dynamics simulations of 74 protein-ligand complexes. Using the calculated binding free energies together with experimental values, we estimated free-energy costs for the large-scale (~17-20 Å) conformational change of the activation loop by an indirect approach, circumventing the very challenging problem of simulating the conformational change directly. We also used the Potts statistical potential to thread large sequence ensembles over active and inactive kinase states. The structure-based and sequence-based analyses are consistent; together they suggest TKs evolved to have free-energy penalties for the classical 'folded activation loop' DFG-out conformation relative to the active conformation, that is, on average, 4-6 kcal/mol smaller than the corresponding values for STKs. Potts statistical energy analysis suggests a molecular basis for this observation, wherein the activation loops of TKs are more weakly 'anchored' against the catalytic loop motif in the active conformation and form more stable substrate-mimicking interactions in the inactive conformation. These results provide insights into the molecular basis for the divergent functional properties of TKs and STKs, and have pharmacological implications for the target selectivity of type-II inhibitors.

AMBER Drug Discovery Boost Tools: Automated Workflow for Production Free-Energy Simulation Setup and Analysis (ProFESSA)

We report an automated workflow for production free-energy simulation setup and analysis (ProFESSA) using the GPU-accelerated AMBER free-energy engine with enhanced sampling features and analysis tools, part of the AMBER Drug Discovery Boost package that has been integrated into the AMBER22 release. The workflow establishes a flexible, end-to-end pipeline for performing alchemical free-energy simulations that brings to bear technologies, including new enhanced sampling features and analysis tools, to practical drug discovery problems. ProFESSA provides the user with top-level control of large sets of free-energy calculations and offers access to the following key functionalities: (1) automated setup of file infrastructure; (2) enhanced conformational and alchemical sampling with the ACES method; and (3) network-wide free-energy analysis with the optional imposition of cycle closure and experimental constraints. The workflow is applied to perform absolute and relative solvation free-energy and relative ligand-protein binding free-energy calculations using different atom-mapping procedures. Results demonstrate that the workflow is internally consistent and highly robust. Further, the application of a new network-wide Lagrange multiplier constraint analysis that imposes key experimental constraints substantially improves binding free-energy predictions.

Forces Driving a Magic Bullet to Its Target: Revisiting the Role of Thermodynamics in Drug Design, Development, and Optimization

Drug discovery strategies have advanced significantly towards prioritizing target selectivity to achieve the longstanding goal of identifying "magic bullets" amongst thousands of chemical molecules screened for therapeutic efficacy. A myriad of emerging and existing health threats, including the SARS-CoV-2 pandemic, alarming increase in bacterial resistance, and potentially fatal chronic ailments, such as cancer, cardiovascular disease, and neurodegeneration, have incentivized the discovery of novel therapeutics in treatment regimens. The design, development, and optimization of lead compounds represent an arduous and time-consuming process that necessitates the assessment of specific criteria and metrics derived via multidisciplinary approaches incorporating functional, structural, and energetic properties. The present review focuses on specific methodologies and technologies aimed at advancing drug development with particular emphasis on the role of thermodynamics in elucidating the underlying forces governing ligand-target interaction selectivity and specificity. In the pursuit of novel therapeutics, isothermal titration calorimetry (ITC) has been utilized extensively over the past two decades to bolster drug discovery efforts, yielding information-rich thermodynamic binding signatures. A wealth of studies recognizes the need for mining thermodynamic databases to critically examine and evaluate prospective drug candidates on the basis of available metrics. The ultimate power and utility of thermodynamics within drug discovery strategies reside in the characterization and comparison of intrinsic binding signatures that facilitate the elucidation of structural-energetic correlations which assist in lead compound identification and optimization to improve overall therapeutic efficacy.

Evaluating the use of absolute binding free energy in the fragment optimisation process

Key to the fragment optimisation process within drug design is the need to accurately capture the changes in affinity that are associated with a given set of chemical modifications. Due to the weakly binding nature of fragments, this has proven to be a challenging task, despite recent advancements in leveraging experimental and computational methods. In this work, we evaluate the use of Absolute Binding Free Energy (ABFE) calculations in guiding fragment optimisation decisions, retrospectively calculating binding free energies for 59 ligands across 4 fragment elaboration campaigns. We first demonstrate that ABFEs can be used to accurately rank fragment-sized binders with an overall Spearman's r of 0.89 and a Kendall τ of 0.67, although often deviating from experiment in absolute free energy values with an RMSE of 2.75 kcal/mol. We then also show that in several cases, retrospective fragment optimisation decisions can be supported by the ABFE calculations. Comparing against cheaper endpoint methods, namely Nwat-MM/GBSA, we find that ABFEs offer better ranking power and correlation metrics. Our results indicate that ABFE calculations can usefully guide fragment elaborations to maximise affinity.

Absolute binding free energy calculations improve enrichment of actives in virtual compound screening

We determined the effectiveness of absolute binding free energy (ABFE) calculations to refine the selection of active compounds in virtual compound screening, a setting where the more commonly used relative binding free energy approach is not readily applicable. To do this, we conducted baseline docking calculations of structurally diverse compounds in the DUD-E database for three targets, BACE1, CDK2 and thrombin, followed by ABFE calculations for compounds with high docking scores. The docking calculations alone achieved solid enrichment of active compounds over decoys. Encouragingly, the ABFE calculations then improved on this baseline. Analysis of the results emphasizes the importance of establishing high quality ligand poses as starting points for ABFE calculations, a nontrivial goal when processing a library of diverse compounds without informative co-crystal structures. Overall, our results suggest that ABFE calculations can play a valuable role in the drug discovery process.

Drug Discovery by Automated Adaptation of Chemical Structure and Identity

Computer-aided drug design offers the potential to dramatically reduce the cost and effort required for drug discovery. While screening-based methods are valuable in the early stages of hit identification, they are frequently succeeded by iterative, hypothesis-driven computations that require recurrent investment of human time and intuition. To increase automation, we introduce a computational method for lead refinement that combines concerted dynamics of the ligand/protein complex via molecular dynamics simulations with integrated Monte Carlo-based changes in the chemical formula of the ligand. This approach, which we refer to as ligand-exchange Monte Carlo molecular dynamics, accounts for solvent- and entropy-based contributions to competitive binding free energies by coupling the energetics of bound and unbound states during the ligand-exchange attempt. Quantitative comparison of relative binding free energies to reference values from free energy perturbation, conducted in vacuum, indicates that ligand-exchange Monte Carlo molecular dynamics simulations sample relevant conformational ensembles and are capable of identifying strongly binding compounds. Additional simulations demonstrate the use of an implicit solvent model. We speculate that the use of chemical graphs in which exchanges are only permitted between ligands with sufficient similarity may enable an automated search to capture some of the benefits provided by human intuition during hypothesis-guided lead refinement.

LIMO: Latent Inceptionism for Targeted Molecule Generation

Generation of drug-like molecules with high binding affinity to target proteins remains a difficult and resource-intensive task in drug discovery. Existing approaches primarily employ reinforcement learning, Markov sampling, or deep generative models guided by Gaussian processes, which can be prohibitively slow when generating molecules with high binding affinity calculated by computationally-expensive physics-based methods. We present Latent Inceptionism on Molecules (LIMO), which significantly accelerates molecule generation with an inceptionism-like technique. LIMO employs a variational autoencoder-generated latent space and property prediction by two neural networks in sequence to enable faster gradient-based reverse-optimization of molecular properties. Comprehensive experiments show that LIMO performs competitively on benchmark tasks and markedly outperforms state-of-the-art techniques on the novel task of generating drug-like compounds with high binding affinity, reaching nanomolar range against two protein targets. We corroborate these docking-based results with more accurate molecular dynamics-based calculations of absolute binding free energy and show that one of our generated drug-like compounds has a predicted K _D (a measure of binding affinity) of 6 · 10^-14 M against the human estrogen receptor, well beyond the affinities of typical early-stage drug candidates and most FDA-approved drugs to their respective targets. Code is available at https://github.com/Rose-STL-Lab/LIMO.

Relative Binding Free Energy between Chemically Distant Compounds Using a Bidirectional Nonequilibrium Approach

In the context of advanced hit-to-lead drug design based on atomistic molecular dynamics simulations, we propose a dual topology alchemical approach for calculating the relative binding free energy (RBFE) between two chemically distant compounds. The method (termed NE-RBFE) relies on the enhanced sampling of the end-states in bulk and in the bound state via Hamiltonian Replica Exchange, alchemically connected by a series of independent and fast nonequilibrium (NE) simulations. The technique has been implemented in a bidirectional fashion, applying the Crooks theorem to the NE work distributions for RBFE predictions. The dissipation of the NE process, negatively affecting accuracy, has been minimized by introducing a smooth regularization based on shifted electrostatic and Lennard-Jones non bonded potentials. As a challenging testbed, we have applied our method to the calculation of the RBFEs in the recent host-guest SAMPL international contest, featuring a macrocyclic host with guests varying in the net charge, volume, and chemical fingerprints. Closure validation has been successfully verified in cycles involving compounds with disparate Tanimoto coefficients, volume, and net charge. NE-RBFE is specifically tailored for massively parallel facilities and can be used with little or no code modification on most of the popular software packages supporting nonequilibrium alchemical simulations, such as Gromacs, Amber, NAMD, or OpenMM. The proposed methodology bypasses most of the entanglements and limitations of the standard single topology RBFE approach for strictly congeneric series based on free-energy perturbation, such as slowly relaxing cavity water, sampling issues along the alchemical stratification, and the need for highly overlapping molecular fingerprints.

Computational Investigations on the Natural Small Molecule as an Inhibitor of Programmed Death Ligand 1 for Cancer Immunotherapy

Several therapeutic monoclonal antibodies approved by the FDA are available against the PD-1/PD-L1 (programmed death 1/programmed death ligand 1) immune checkpoint axis, which has been an unprecedented success in cancer treatment. However, existing therapeutics against PD-L1, including small molecule inhibitors, have certain drawbacks such as high cost and drug resistance that challenge the currently available anti-PD-L1 therapy. Therefore, this study presents the screening of 32,552 compounds from the Natural Product Atlas database against PD-L1, including three steps of structure-based virtual screening followed by binding free energy to refine the ideal conformation of potent PD-L1 inhibitors. Subsequently, five natural compounds, i.e., Neoenactin B1, Actinofuranone I, Cosmosporin, Ganocapenoid A, and 3-[3-hydroxy-4-(3-methylbut-2-enyl)phenyl]-5-(4-hydroxybenzyl)-4-methyldihydrofuran-2(3H)-one, were collected based on the ADMET (absorption, distribution, metabolism, excretion, and toxicity) profiling and binding free energy (>−60 kcal/mol) for further computational investigation in comparison to co-crystallized ligand, i.e., JQT inhibitor. Based on interaction mapping, explicit 100 ns molecular dynamics simulation, and end-point binding free energy calculations, the selected natural compounds were marked for substantial stability with PD-L1 via intermolecular interactions (hydrogen and hydrophobic) with essential residues in comparison to the JQT inhibitor. Collectively, the calculated results advocate the selected natural compounds as the putative potent inhibitors of PD-L1 and, therefore, can be considered for further development of PD-L1 immune checkpoint inhibitors in cancer immunotherapy.

Gene Interaction Network Analysis Reveals IFI44L as a Drug Target in Rheumatoid Arthritis and Periodontitis

Simple Summary

In spite of substantial investigation, the biological link between periodontitis and rheumatoid arthritis remains unexplained. This study intended to correlate periodontitis and rheumatoid arthritis gene expression patterns to find shared targets for both the disease. We identified the differentially expressed genes (DEGs) in periodontitis and rheumatoid arthritis. The network was built by integrating DEGs and ranking the genes using GeneMANIA. FINDSITEcomb2.0 was used to find a possible inhibitor for the top-ranked gene. Further, the binding effectiveness and protein-ligand complex stability were then determined by molecular docking and molecular dynamics. The network analysis showed IFI44L as a highly ranking gene implicated in most immunological pathways. A virtual screening of 6507 compounds revealed vemurafenib as the best candidate for the IFI44L target. Molecular docking and molecular dynamics modelling revealed the stability of the IFI44L-vemurafenib complex, which suggest IFI44L is potential target and vemurafenib could be the better candidate to treat both diseases.

Abstract

Objective: Despite extensive research on periodontitis and rheumatoid arthritis, the underlying molecular connectivity between these condition remains largely unknown. This research aimed to integrate periodontitis and rheumatoid arthritis gene expression profiles to identify interconnecting genes and focus to develop a common lead molecule against these inflammatory conditions. Materials and Methods: Differentially expressed genes (DEGs) of periodontitis and rheumatoid arthritis were identified from the datasets retrieved from the Gene Expression Omnibus database. The network was constructed by merging DEGs, and the interconnecting genes were identified and ranked using GeneMANIA. For the selected top ranked gene, the potential inhibitor was searched using FINDSITE^comb2.0. Subsequently, the molecular docking and molecular dynamics were performed to determine the binding efficiency and protein-ligand complex stability, respectively. Results: From the network analysis, IFN-induced protein 44-like (IFI44L) was identified as a top ranked gene involved in most of the immunological pathway. With further virtual screening of 6507 molecules, vemurafenib was identified to be the best fit against the IFI44L target. The binding energy and stability of IFI44L with vemurafenib were investigated using molecular docking and molecular dynamics simulation. Docking results show binding energy of −7.7 Kcal/mol, and the simulation results show stability till 100 ns. Conclusions: The identified IFI44L may represent a common drug target for periodontitis and rheumatoid arthritis. Vemurafenib could be a potent anti-inflammatory drug for both diseases.

Potential Therapeutic Targets of Rehmannia Formulations on Diabetic Nephropathy: A Comparative Network Pharmacology Analysis

Background: Previous retrospective cohorts showed that Rehmannia-6 (R-6, Liu-wei-di-huang-wan) formulations were associated with significant kidney function preservation and mortality reduction among chronic kidney disease patients with diabetes. This study aimed to investigate the potential mechanism of action of common R-6 variations in a clinical protocol for diabetic nephropathy (DN) from a system pharmacology approach.

Study Design and Methods: Disease-related genes were retrieved from GeneCards and OMIM by searching “Diabetic Nephropathy” and “Macroalbuminuria”. Variations of R-6 were identified from a published existing clinical practice guideline developed from expert consensus and pilot clinical service program. The chemical compound IDs of each herb were retrieved from TCM-Mesh and PubChem. Drug targets were subsequently revealed via PharmaMapper and UniProtKB. The disease gene interactions were assessed through STRING, and disease–drug protein–protein interaction network was integrated and visualized by Cytoscape. Clusters of disease–drug protein–protein interaction were constructed by Molecular Complex Detection (MCODE) extension. Functional annotation of clusters was analyzed by DAVID and KEGG pathway enrichment. Differences among variations of R-6 were compared. Binding was verified by molecular docking with AutoDock.

Results: Three hundred fifty-eight genes related to DN were identified, forming 11 clusters which corresponded to complement and coagulation cascades and signaling pathways of adipocytokine, TNF, HIF-1, and AMPK. Five variations of R-6 were analyzed. Common putative targets of the R-6 variations on DN included ACE, APOE, CCL2, CRP, EDN1, FN1, HGF, ICAM1, IL10, IL1B, IL6, INS, LEP, MMP9, PTGS2, SERPINE1, and TNF, which are related to regulation of nitric oxide biosynthesis, lipid storage, cellular response to lipopolysaccharide, inflammatory response, NF-kappa B transcription factor activity, smooth muscle cell proliferation, blood pressure, cellular response to interleukin-1, angiogenesis, cell proliferation, peptidyl-tyrosine phosphorylation, and protein kinase B signaling. TNF was identified as the seed for the most significant cluster of all R-6 variations. Targets specific to each formulation were identified. The key chemical compounds of R-6 have good binding ability to the putative protein targets.

Conclusion: The mechanism of action of R-6 on DN is mostly related to the TNF signaling pathway as a core mechanism, involving amelioration of angiogenesis, fibrosis, inflammation, disease susceptibility, and oxidative stress. The putative targets identified could be validated through clinical trials.

Comparing the Performance of Different AMBER Protein Forcefields, Partial Charge Assignments, and Water Models for Absolute Binding Free Energy Calculations

Identifying chemical starting points is a vital first step in small molecule drug discovery and can take significant time and money. For this reason, computational approaches to virtual screening are of great interest as they can lower the cost and shorten timeframes. However, simple approaches such as molecular docking and pharmacophore screening are of limited accuracy and provide a low probability of success. Alchemical binding free energies represent a promising approach for virtual screening as they naturally incorporate the key effects of water molecules, protein flexibility, and binding entropy. However, the calculations are technically very challenging, with performance depending on the specific forcefield used. For this reason, it is important that the community has access to benchmark test sets to assess prediction accuracy. In this paper, we present an approach to alchemical binding free energies using OpenMM. We identify effective simulation parameters using an existing BRD4(1) test set and present two new benchmark sets (cMET and PDE2A) that can be used in the community for validation purposes. Our findings also highlight the effectiveness of some AMBER forcefields, in particular, AMBER ff15ipq.

Structure-based virtual screening workflow to identify antivirals targeting HIV-1 capsid

We have identified novel HIV-1 capsid inhibitors targeting the PF74 binding site. Acting as the building block of the HIV-1 capsid core, the HIV-1 capsid protein plays an important role in the viral life cycle and is an attractive target for antiviral development. A structure-based virtual screening workflow for hit identification was employed, which includes docking 1.6 million commercially-available drug-like compounds from the ZINC database to the capsid dimer, followed by applying two absolute binding free energy (ABFE) filters on the 500 top-ranked molecules from docking. The first employs the Binding Energy Distribution Analysis Method (BEDAM) in implicit solvent. The top-ranked compounds are then refined using the Double Decoupling method in explicit solvent. Both docking and BEDAM refinement were carried out on the IBM World Community Grid as part of the FightAIDS@Home project. Using this virtual screening workflow, we identified 24 molecules with calculated binding free energies between − 6 and − 12 kcal/mol. We performed thermal shift assays on these molecules to examine their potential effects on the stability of HIV-1 capsid hexamer and found that two compounds, ZINC520357473 and ZINC4119064 increased the melting point of the latter by 14.8 °C and 33 °C, respectively. These results support the conclusion that the two ZINC compounds are primary hits targeting the capsid dimer interface. Our simulations also suggest that the two hit molecules may bind at the capsid dimer interface by occupying a new sub-pocket that has not been exploited by existing CA inhibitors. The possible causes for why other top-scored compounds suggested by ABFE filters failed to show measurable activity are discussed.

Developing end-point methods for absolute binding free energy calculation using the Boltzmann-quasiharmonic model

Understanding the physical forces underlying receptor-ligand binding requires robust methods for analyzing the binding thermodynamics. In end-point binding free energy methods the binding free energy is naturally decomposable into physically intuitive contributions such as the solvation free energy and configurational entropy that can provide insights. Here we present a new end-point method called EE-BQH (Effective Energy-Boltzmann-Quasiharmonic) which combines the Boltzmann-Quasiharmonic model for configurational entropy with different solvation free energy methods, such as the continuum solvent PBSA model and the integral equation-based 3D-RISM, to estimate the absolute binding free energy. We compare EE-BQH with other treatments of configurational entropy such as Quasiharmonic models in internal coordinates (QHIC) and in Cartesian coordinates (QHCC), and Normal Mode analysis (NMA), by testing them on the octa acids host-guest complexes from the SAMPL8 blind challenge. The accuracies in the calculated absolute binding free energies strongly depend on the configurational entropy and solvation free energy methods used. QHIC and BQH yield the best agreements with the established potential of mean force (PMF) estimates, with R² of ∼0.7 and mean unsigned error of ∼1.7 kcal mol^-1. These results from the end-point calculations are also in similar agreement with experiments. While 3D-RISM in combination with QHIC or BQH lead to reasonable correlations with the PMF results and experiments, the calculated absolute binding free energies are underestimated by ∼5 kcal mol^-1. While the binding is accompanied by a significant reduction in the ligand translational/rotational entropy, the change in the torsional entropy in these host-guest systems is slightly positive. Compared with BQH, QHIC underestimates the reduction of configurational entropy because of the non-Gaussian probability distributions in the ligand rotation and a small number of torsions. The study highlights the crucial role of configurational entropy in determining binding and demonstrates the potential of using the new end-point method to provide insights in more complex protein-ligand systems.

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

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

Symmetry-Adapted Restraints for Binding Free Energy Calculations

Binding free energy calculations rely critically on a precise definition of the bound state and well-designed ligand restraints to ensure that binding free energy calculations converge rapidly and yield estimates of well-defined thermodynamic quantities. The distance-to-bound-configuration (DBC) is a single variable that can precisely delineate the bound state of a ligand including translational, rotational and conformational degrees of freedom and has been successfully used to capture binding modes with complex geometries. DBC is defined as the root-mean-square deviation (RMSD) of ligand coordinates in the frame of reference of the binding site. In the special case where the ligand features symmetry-equivalent atoms, a standard RMSD arbitrarily distinguishes equivalent poses, mixing equivalent and nonequivalent degrees of freedom, and preventing the precise delineation of the bound state ensemble, which negates the benefits of defining a flat-bottom binding restraint. To remedy this, we introduce a symmetry-adapted DBC coordinate where the RMSD is minimized over permutations of equivalent ligand atoms. This coordinate is implemented in a portable software library, the Collective Variables Module. We tested the approach by computing the absolute binding free energy of benzene to the engineered site of a mutant lysozyme (L99A/M102H) using alchemical free energy perturbation. We found that the symmetry-adapted restraint leads to well-behaved convergence of both the decoupling free energy in the binding site and the restrained free energy in the gas phase, recovering the affinity computed using a classic center-of-mass restraint. Thus, symmetry-adapted DBC seamlessly generalizes the benefits of DBC restraints to the case of symmetric ligands. The underlying symmetric RMSD coordinate can also be used for analyzing or biasing simulations in other contexts than affinity predictions.

Fast Equilibration of Water between Buried Sites and the Bulk by Molecular Dynamics with Parallel Monte Carlo Water Moves on Graphical Processing Units

Molecular dynamics (MD) simulations of proteins are commonly used to sample from the Boltzmann distribution of conformational states, with wide-ranging applications spanning chemistry, biophysics, and drug discovery. However, MD can be inefficient at equilibrating water occupancy for buried cavities in proteins that are inaccessible to the surrounding solvent. Indeed, the time needed for water molecules to equilibrate between the bulk solvent and the binding site can be well beyond what is practical with standard MD, which typically ranges from hundreds of nanoseconds to a few microseconds. We recently introduced a hybrid Monte Carlo/MD (MC/MD) method, which speeds up the equilibration of water between buried cavities and the surrounding solvent, while sampling from the thermodynamically correct distribution of states. While the initial implementation of the MC functionality led to considerable slowing of the overall simulations, here we address this problem with a parallel MC algorithm implemented on graphical processing units. This results in speed-ups of 10-fold to 1000-fold over the original MC/MD algorithm, depending on the system and simulation parameters. The present method is available for use in the AMBER simulation software.