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

Lambda-ABF-OPES: Faster Convergence with High Accuracy in Alchemical Free Energy Calculations

Predicting the binding affinity between small molecules and target macromolecules while combining both speed and accuracy is a cornerstone of modern computational drug discovery, which is critical for accelerating therapeutic development. Despite recent progress in molecular dynamics (MD) simulations, such as advanced polarizable force fields and enhanced sampling techniques, estimating absolute binding free energies (ABFEs) remains computationally challenging. To overcome these difficulties, we introduce a highly efficient hybrid methodology that couples the Lambda-adaptive biasing force (Lambda-ABF) scheme with on-the-fly probability enhanced sampling (OPES). This approach achieves up to a 9-fold improvement in sampling efficiency and computational speed compared to the original Lambda-ABF when used in conjunction with the AMOEBA polarizable force field, yielding converged results at a fraction of the cost of standard techniques.

Computation of Protein-Ligand Binding Free Energies with a Quantum Mechanics-Based Mining Minima Algorithm

A new method, protein-ligand QM-VM2 (PLQM-VM2), to calculate protein-ligand binding free energies by combining mining minima, a statistical mechanics end-point-based approach, with quantum mechanical potentials is presented. PLQM-VM2 is described in terms of a highly flexible workflow that is initiated from a Protein Data Bank (PDB) file and a chemical structure data file (SD file) containing two-dimensional (2D) or three-dimensional (3D) ligand series coordinates. The workflow utilizes the previously developed molecular mechanics (MM) implementation of the second-generation mining minima method, MM-VM2, to provide ensembles of protein, free ligand, and protein-ligand conformers, which are postprocessed at chosen levels of QM theory, via the quantum chemistry software package GAMESS, to correct MM-based conformer geometries and electronic energies. The corrected energies are used in the calculation of configuration integrals, which on summation over the conformer ensembles give QM-corrected chemical potentials and ultimately QM-corrected binding free energies. In this work, PLQM-VM2 is applied to three benchmark protein-ligand series: HIV-1 protease/38 ligands, c-Met/24 ligands, and TNKS2/27 ligands. QM corrections are carried out at the semiempirical third-order density functional tight-binding level of theory, augmented with dispersion and damping corrections (DFTB3-D3(BJ)H). Bulk solvation effects are accounted for with the conductor-like polarizable continuum model (PCM). DFTB3-D3(BJ)H/PCM single-point energy-only and geometry optimization QM corrections are carried out in conjunction with two different models that address the large computational scaling associated with protein-sized molecular systems. One is a protein cutout model, whereby a set of protein atoms in and around the binding site are carved out, dangling bonds are capped with hydrogens, and only atoms directly in the protein binding site are mobile along with the ligand atoms. The other model is the Fragment Molecular Orbital (FMO) method, which includes the whole protein system but again only allows the binding site and ligand atoms to be mobile. All four of these methodological approaches to QM corrections provide significant improvement over MM-VM2 in terms of rank order and parametric linear correlation with experimentally determined binding affinities. Overall, FMO with geometry optimizations performed the best, but the much cheaper cutout single-point energy approach still provides a good level of accuracy. Furthermore, a clear result is that the PLQM-VM2 calculated binding free energies for the three diverse test systems in this work are, in contrast to those calculated using MM-VM2, directly comparable in energy scale. This suggests a basis for future development of a PLQM-VM2-based multiprotein screening capability to check for off-target activity of ligand series. Benchmark timings on a single compute node (32 CPU cores) for PLQM-VM2 calculation of the chemical potential of a protein-ligand complex range from ca. 30-45 min for the single-point energy approaches to ∼5 h for the cutout approach with geometry optimization and to ∼35 h for the full protein FMO approach with geometry optimization.

Relative Binding Free Energy Estimation of Congeneric Ligands and Macromolecular Mutants with the Alchemical Transfer Method with Coordinate Swapping

We present the Alchemical Transfer with Coordinate Swapping (ATS) method to enable the calculation of the relative binding free energies between large congeneric ligands and single-point mutant peptides to protein receptors with the Alchemical Transfer Method (ATM) framework. Similarly to ATM, the new method implements the alchemical transformation as a coordinate transformation and works with any unmodified force fields and standard chemical topologies. Unlike ATM, which transfers whole ligands in and out of the receptor binding site, ATS limits the magnitude of the alchemical perturbation by transferring only the portion of the molecules that differ between the bound and unbound ligands. The common region of the two ligands, which can be arbitrarily large, is unchanged and does not contribute to the magnitude and statistical fluctuations of the perturbation energy. Internally, the coordinates of the atoms of the common regions are swapped to maintain the integrity of the covalent bonding data structures of the OpenMM molecular dynamics engine. The work successfully validates the method on protein-ligand and protein-peptide RBFE benchmarks. This advance paves the road for the application of the relative binding free energy Alchemical Transfer Method protocol to study the effect of protein and nucleic acid mutations on the binding affinity and specificity of macromolecular complexes.

Predicting Thermodynamic Stability at Protein G Sites with Deleterious Mutations Using λ-Dynamics with Competitive Screening

Free energy predictions are useful in protein design and computer-aided drug design. Alchemical free energy methods are highly accurate, and the alchemical method λ-dynamics significantly improves computational cost. Recent progress made simulations of dozens of perturbations at a single site possible, enabling in silico site-saturation mutagenesis with λ-dynamics. Site-saturation mutagenesis may require increased sampling to characterize many mutations and to accommodate structural disruptions around deleterious mutations. We reintroduce the neglected idea of competitive screening with λ-dynamics to address both issues. Traditional landscape flattening tunes two distinct biases to sample all mutations equally in the folded and unfolded states. Competitive screening transfers the unfolded bias to the folded state so that only reasonable mutations are sampled. Competitive screening is demonstrated on four surface sites and four buried sites in protein G and provides improvements for buried sites. Consequently, competitive screening provides new opportunities for molecular design within larger chemical spaces.

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.

Robust protein-ligand interaction modeling through integrating physical laws and geometric knowledge for absolute binding free energy calculation

Accurate estimation of protein-ligand (PL) binding free energies is a crucial task in medicinal chemistry and a critical measure of PL interaction modeling effectiveness. However, traditional computational methods are often computationally expensive and prone to errors. Recently, deep learning (DL)-based approaches for predicting PL interactions have gained enormous attention, but their accuracy and generalizability are hindered by data scarcity. In this study, we propose LumiNet, a versatile PL interaction modeling framework that bridges the gap between physics-based models and black-box algorithms. LumiNet utilizes a subgraph transformer to extract multiscale information from molecular graphs and employs geometric neural networks to integrate PL information, mapping atomic pair structures into key physical parameters of non-bonded interactions in classical force fields, thereby enhancing accurate absolute binding free energy (ABFE) calculations. LumiNet is designed to be highly interpretable, offering detailed insights into atomic interactions within protein-ligand complexes, pinpointing relatively important atom pairs or groups. Our semi-supervised learning strategy enables LumiNet to adapt to new targets with fewer data points than other data-driven methods, making it more relevant for real-world drug discovery. Benchmarks show that LumiNet outperforms the current state-of-the-art model by 18.5% on the PDE10A dataset, and rivals the FEP+ method in some tests with a speed improvement of several orders of magnitude. We applied LumiNet in the scaffold hopping process, which accurately guided the discovery of the optimal ligands. Furthermore, we provide a web service for the research community to test LumiNet. The visualization of predicted inter-molecular energy contributions is expected to provide practical value in drug discovery projects.

The Need for Continuing Blinded Pose- and Activity Prediction Benchmarks

Computational tools for structure-based drug design (SBDD) are widely used in drug discovery and can provide valuable insights to advance projects in an efficient and cost-effective manner. However, despite the importance of SBDD to the field, the underlying methodologies and techniques have many limitations. In particular, binding pose and activity predictions (P-AP) are still not consistently reliable. We strongly believe that a limiting factor is the lack of a widely accepted and established community benchmarking process that independently assesses the performance and drives the development of methods, similar to the CASP benchmarking challenge for protein structure prediction. Here, we provide an overview of P-AP, unblinded benchmarking data sets, and blinded benchmarking initiatives (concluded and ongoing) and offer a perspective on learnings and the future of the field. To accelerate a breakthrough on the development of novel P-AP methods, it is necessary for the community to establish and support a long-term benchmark challenge that provides nonbiased training/test/validation sets, a systematic independent validation, and a forum for scientific discussions.

IMERGE-FEP: Improving Relative Free Energy Calculation Convergence with Chemical Intermediates

Alchemical free energy calculations are becoming an increasingly prevalent tool in drug discovery efforts. Over the past decade, significant progress has been made in automating various aspects of this technique. However, one aspect hampering wider application is the construction of perturbation networks to connect ligands of interest. More specifically, ligand pairs with large dissimilarities should be avoided since they can lower convergence and decrease accuracy. Here, we propose a technique for automatic generation of intermediate molecules to break up problematic edges─calculations connecting two different ligands or molecules─into smaller perturbations. To this end, a modular tool was developed that generates intermediates for a molecule pair by enumerating R-group combinations called IMERGE-FEP (Intermediate MolEculaR GEnerator for Free Energy Perturbation). Intermediate enumeration of multiple, representative congeneric series showed that intermediates increase similarity regarding shared substructures, geometry, and LOMAP scores. Taken together, this tool eases integration of intermediate steps into free energy calculation protocols.

High performance-oriented computer aided drug design approaches in the exascale era

INTRODUCTION

In 2023, the first exascale supercomputer was opened to the public in the US. With a demonstrated 1.1 exaflops of performance, Frontier represents an unprecedented breakthrough in high-performance computing (HPC). Currently, more (and more powerful) machines are being installed worldwide. Computer-aided drug design (CADD) is one of the fields of computational science that can greatly benefit from exascale computing for the benefit of the whole society. However, scaling CADD approaches to exploit exascale machines require new algorithmic and software solutions.

AREAS COVERED

Here, the authors consider physics-based and machine learning (ML)-aided techniques for the design of small molecule binders capable of leveraging modern parallel computer architectures. Specifically, the authors focus on HPC-oriented large-scale applications from the past 3 years that were enabled by (pre)exascale supercomputers by running on up tothousands of accelerated nodes.

EXPERT OPINION

In the area of ML, exascale computers can enable the training of generative models with unprecedented predictive power to design novel ligands, provided large amounts of high-quality data are available. Exascale computers could also unlock the potential of accurate ML-aided physics-based methods to boost the success rate of structure-based drug design campaigns. Currently, however, methodological developments are still required to allow routine large-scale applications of such rigorous approaches.

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.

Normalized Protein-Ligand Distance Likelihood Score for End-to-End Blind Docking and Virtual Screening

Molecular Docking is a critical task in structure-based virtual screening. Recent advancements have showcased the efficacy of diffusion-based generative models for blind docking tasks. However, these models do not inherently estimate protein-ligand binding strength thus cannot be directly applied to virtual screening tasks. Protein-ligand scoring functions serve as fast and approximate computational methods to evaluate the binding strength between the protein and ligand. In this work, we introduce normalized mixture density network (NMDN) score, a deep learning (DL)-based scoring function learning the probability density distribution of distances between protein residues and ligand atoms. The NMDN score addresses limitations observed in existing DL scoring functions and performs robustly in both pose selection and virtual screening tasks. Additionally, we incorporate an interaction module to predict the experimental binding affinity score to fully utilize the learned protein and ligand representations. Finally, we present an end-to-end blind docking and virtual screening protocol named DiffDock-NMDN. For each protein-ligand pair, we employ DiffDock to sample multiple poses, followed by utilizing the NMDN score to select the optimal binding pose, and estimating the binding affinity using scoring functions. Our protocol achieves an average enrichment factor of 4.96 on the LIT-PCBA data set, proving effective in real-world drug discovery scenarios where binder information is limited. This work not only presents a robust DL-based scoring function with superior pose selection and virtual screening capabilities but also offers a blind docking protocol and benchmarks to guide future scoring function development.

Narrowing the gap between machine learning scoring functions and free energy perturbation using augmented data

Machine learning offers great promise for fast and accurate binding affinity predictions. However, current models lack robust evaluation and fail on tasks encountered in (hit-to-) lead optimisation, such as ranking the binding affinity of a congeneric series of ligands, thereby limiting their application in drug discovery. Here, we address these issues by first introducing a novel attention-based graph neural network model called AEV-PLIG (atomic environment vector-protein ligand interaction graph). Second, we introduce a new and more realistic out-of-distribution test set called the OOD Test. We benchmark our model on this set, CASF-2016, and a test set used for free energy perturbation (FEP) calculations, that not only highlights the competitive performance of AEV-PLIG, but provides a realistic assessment of machine learning models with rigorous physics-based approaches. Moreover, we demonstrate how leveraging augmented data (generated using template-based modelling or molecular docking) can significantly improve binding affinity prediction correlation and ranking on the FEP benchmark (weighted mean PCC and Kendall's τ increases from 0.41 and 0.26 to 0.59 and 0.42). These strategies together are closing the performance gap with FEP calculations (FEP+ achieves weighted mean PCC and Kendall's τ of 0.68 and 0.49 on the FEP benchmark) while being  ~400,000 times faster.

Potential distribution theory of alchemical transfer

We present an analytical description of the Alchemical Transfer Method (ATM) for molecular binding using the Potential Distribution Theory (PDT) formalism. ATM models the binding free energy by mapping the bound and unbound states of the complex by translating the ligand coordinates. PDT relates the free energy and the probability densities of the perturbation energy along the alchemical path to the probability density at the initial state, which is the unbound state of the complex in the case of a binding process. Hence, the ATM probability density of the transfer energy at the unbound state is first related by a convolution operation of the probability densities for coupling the ligand to the solvent and coupling it to the solvated receptor-for which analytical descriptions are available-with parameters obtained from maximum likelihood analysis of data from double-decoupling alchemical calculations. PDT is then used to extend this analytical description along the alchemical transfer pathway. We tested the theory on the alchemical binding of five guests to the tetramethyl octa-acid host from the SAMPL8 benchmark set. In each case, the probability densities of the perturbation energy for transfer along the alchemical transfer pathway obtained from numerical calculations match those predicted from the theory and double-decoupling simulations. The work provides a solid theoretical foundation for alchemical transfer, offers physical insights on the form of the probability densities observed in alchemical transfer calculations, and confirms the conceptual and numerical equivalence between the alchemical transfer and double-decoupling processes.

Evaluation of Machine Learning/Molecular Mechanics End-State Corrections with Mechanical Embedding to Calculate Relative Protein-Ligand Binding Free Energies

The development of machine-learning (ML) potentials offers significant accuracy improvements compared to molecular mechanics (MM) because of the inclusion of quantum-mechanical effects in molecular interactions. However, ML simulations are several times more computationally demanding than MM simulations, so there is a trade-off between speed and accuracy. One possible compromise are hybrid machine learning/molecular mechanics (ML/MM) approaches with mechanical embedding that treat the intramolecular interactions of the ligand at the ML level and the protein-ligand interactions at the MM level. Recent studies have reported improved protein-ligand binding free energy results based on ML/MM using ANI-2x with mechanical embedding, arguing that intramolecular interactions like torsion potentials of the ligand are often the limiting factor for accuracy. This claim is evaluated based on 108 relative binding free energy calculations for four different benchmark systems. As an alternative strategy, we also tested a tool that fits the MM dihedral potentials to the ML level of theory. Fitting was performed with the ML potentials ANI-2x and AIMNet2, and, for the benchmark system TYK2, also with quantum-mechanical calculations using ωB97M-D3(BJ)/def2-TZVPPD. Overall, the relative binding free energy results from MM with Open Force Field 2.2.0, MM with ML-fitted torsion potentials, and the corresponding ML/MM end-state corrected simulations show no statistically significant differences in the mean absolute errors (between 0.8 and 0.9 kcal mol-1). This can probably be explained by the usage of the same MM parameters to calculate the protein-ligand interactions. Therefore, a well-parametrized force field is on a par with simple mechanical embedding ML/MM simulations for protein-ligand binding. In terms of computational costs, the reparametrization of poor torsional potentials is preferable over employing computationally intensive ML/MM simulations of protein-ligand complexes with mechanical embedding. Also, the refitting strategy leads to lower variances of the protein-ligand binding free energy results than the ML/MM end-state corrections. For free energy corrections with ML/MM, the results indicate that better convergence and more advanced ML/MM schemes will be required for applications in computer-guided drug discovery.

PairMap: An Intermediate Insertion Approach for Improving the Accuracy of Relative Free Energy Perturbation Calculations for Distant Compound Transformations

Accurate prediction of the difference in binding free energy between compounds is crucial for reducing the high costs associated with drug discovery. Relative binding free energy perturbation (RBFEP) calculations are effective for small structural changes; however, large topological changes pose significant challenges for calculations, leading to high errors and difficulties in convergence. To address such issues, we propose a new approach─PairMap─that focuses on introducing appropriate intermediates for complex transformations between two input compounds. PairMap-generated intermediates exhaustively, determined the optimal conversion paths, and introduced thermodynamic cycles into the perturbation map to improve accuracy and reduce computational cost. PairMap succeeded in introducing appropriate intermediates that could not be discovered by existing simple approaches by comprehensively considering intermediates. Furthermore, we evaluated the accuracy of the prediction of binding free energy using 9 compounds selected from Wang et al.'s benchmark set, which included particularly complex transformations. The perturbation map generated by PairMap achieved excellent accuracy with a mean absolute error of 0.93 kcal/mol compared to 1.70 kcal/mol when using the perturbation map generated by the conventional Flare FEP intermediate introduction method. Moreover, in a scaffold hopping experiment conducted with the PDE5a target involving complex transformations, PairMap provided more accurate free energy predictions than ABFEP calculations, yielding more reliable results compared to experimental data. Additionally, PairMap can be utilized to introduce intermediates into congeneric series, demonstrating that complex links on the perturbation map can be resolved with minimal addition of intermediates and links. In conclusion, PairMap overcomes the limitations of existing methods by enabling RBFEP calculations for more complex transformations, further streamlining lead optimization in drug discovery.

Leveraging the Thermodynamics of Protein Conformations in Drug Discovery

As the name implies, structure-based drug design requires confidence in the holo complex structure. The ability to clarify which protein conformation to use when ambiguity arises would be incredibly useful. We present a large scale validation of the computational method Protein Reorganization Free Energy Perturbation (PReorg-FEP) and demonstrate its quantitative accuracy in selecting the correct protein conformation among candidate models in apo or ligand induced states for 14 different systems. These candidate conformations are pulled from various drug discovery related campaigns: cryptic conformations induced by novel hits in lead identification, binding site rearrangement during lead optimization, and conflicting structural biology models. We also show an example of a pH-dependent conformational change, relevant to protein design.

RED-E-Function-Based Equilibrium Parameter Finder: Finding the Best Restraint Parameters in Absolute Binding Free Energy Calculations

Free energy perturbation (FEP)-based absolute binding free energy (ABFE) calculations have emerged as a powerful tool for the accurate prediction of ligand-protein binding affinities in drug discovery. The restraint addition is crucial in FEP-ABFE calculations; however, due to the non-orthogonal couplings between the restrained degrees of freedom, it typically requires numerous λ windows to ensure the phase-space overlap during restraint addition. This study introduces the RED-E-function-based equilibrium parameter finder (REPF), a novel method that relies on harmonic restraints to optimize the equilibrium values in restraints, enhancing phase-space overlap and improving the convergence of the restraint addition. REPF was applied to 44 protein-ligand complexes across 5 targets and compared to restraint schemes reported in the literature. We found that REPF-optimized restraints achieve an accuracy comparable to that of the 12λ approach while using only 2λ simulations, resulting in a significant reduction in computational costs. Extensive tests confirmed the improved convergence behavior and reduced energy fluctuations of REPF-optimized restraints.

Fine-tuning molecular mechanics force fields to experimental free energy measurements

Alchemical free energy methods using molecular mechanics (MM) force fields are essential tools for predicting thermodynamic properties of small molecules, especially via free energy calculations that can estimate quantities relevant for drug discovery such as affinities, selectivities, the impact of target mutations, and ADMET properties. While traditional MM forcefields rely on hand-crafted, discrete atom types and parameters, modern approaches based on graph neural networks (GNNs) learn continuous embedding vectors that represent chemical environments from which MM parameters can be generated. Excitingly, GNN parameterization approaches provide a fully end-to-end differentiable model that offers the possibility of systematically improving these models using experimental data. In this study, we treat a pretrained GNN force field-here, espaloma-0.3.2-as a foundation simulation model and fine-tune its charge model using limited quantities of experimental hydration free energy data, with the goal of assessing the degree to which this can systematically improve the prediction of other related free energies. We demonstrate that a highly efficient "one-shot fine-tuning" method using an exponential (Zwanzig) reweighting free energy estimator can improve prediction accuracy without the need to resimulate molecular configurations. To achieve this "one-shot" improvement, we demonstrate the importance of using effective sample size (ESS) regularization strategies to retain good overlap between initial and fine-tuned force fields. Moreover, we show that leveraging low-rank projections of embedding vectors can achieve comparable accuracy improvements as higher-dimensional approaches in a variety of data-size regimes. Our results demonstrate that linearly-perturbative fine-tuning of foundation model electrostatic parameters to limited experimental data offers a cost-effective strategy that achieves state-of-the-art performance in predicting hydration free energies on the FreeSolv dataset.

Fuentes L, Pimenta A, Nagem RAP, Chávez-Olórtegui C, Schneider FS, Molina F, Sanchez EF, Suárez D, Ferreira RS. Assessing the Interactions between Snake Venom Metalloproteinases and Hydroxamate Inhibitors Using Kinetic and ITC Assays, Molecular Dynamics Simulations and MM/PBSA-Based Scoring Functions

Bothrops species are the main cause of snake bites in rural communities of tropical developing countries of Central and South America. Envenomation by Bothrops snakes is characterized by prominent local inflammation, hemorrhage and necrosis as well as systemic hemostatic disturbances. These pathological effects are mainly caused by the major toxins of the viperidae venoms, the snake venom metalloproteinases (SVMPs). Despite the antivenom therapy efficiency to block the main toxic effects on bite victims, this treatment shows limited efficacy to prevent tissue necrosis. Thus, drug-like inhibitors of these toxins have the potential to aid serum therapy of accidents inflicted by viper snakes. Broad-spectrum metalloprotease inhibitors bearing a hydroxamate zinc-binding group are potential candidates to improve snake bites therapy and could also be used to study toxin-ligand interactions. Therefore, in this work, we used both docking calculations and molecular dynamics simulations to assess the interactions between six hydroxamate inhibitors and two P-I SVMPs selected as models: Atroxlysin-I (hemorrhagic) from Bothrops atrox, and Leucurolysin-a (nonhemorrhagic) from Bothrops leucurus. We also employed a large variety of end-point free energy methods in combination with entropic terms to produce scoring functions of the relative affinities of the inhibitors for the toxins. Then we identified the scoring functions that best correlated with experimental data obtained from kinetic activity assays. In addition, to the characterization of these six molecules as inhibitors of the toxins, this study sheds light on the main enzyme-inhibitor interactions, explaining the broad-spectrum behavior of the inhibitors, and identifies the energetic and entropic terms that improve the performance of the scoring functions.

Protocols for Metallo- and Serine-β-Lactamase Free Energy Predictions: Insights from Cross-Class Inhibitors

While relative binding free energy (RBFE) calculations using alchemical methods are routinely carried out for many pharmaceutically relevant protein targets, challenges remain. For example, open-source tools do not support the easy setup and simulation of metalloproteins, particularly when ligands directly coordinate to the metal site. Here, we evaluate the performance of RBFE methods for KPC-2, a serine-β-lactamase (SBL), and two nonbonded metal parameter setups for VIM-2, a metallo-β-lactamase (MBL) with two active site zinc ions. We tested two different ways of modeling the ligand-zinc interactions. First, a restraint-based approach, in which FF14SB zinc parameters are combined with harmonic restraints between the zincs and their coordinating residues. The second approach uses an upgraded Amber force field (UAFF) for zinc-metalloproteins with adjusted partial charges and nonbonded terms of zinc-coordinating residues. Molecular mechanics (MM) and quantum mechanics/molecular mechanics (QM/MM) simulations show that the crystallographically observed zinc coordination is not retained in MM simulations with either zinc parameter set for a series of known phosphonic acid-based inhibitors bound to VIM-2. These phosphonic acid-based inhibitors exhibit known cross-class affinity for SBLs and MBLs and serve as a benchmark for RBFE calculations for VIM-2, after validation with KPC-2. The KPC-2 free energy of binding estimates are within expected literature accuracies for the ligand series with a mean absolute error of 0.45 0.28 0.66 kcal/mol and a Pearson's correlation coefficient of 0.93 0.85 0.98 . For VIM-2, the UAFF approach has improved correlation from 0.55 - 0.04 0.88 to 0.78 0.38 0.92 , compared to the restraint approach. The presented strategies for handling ligands coordinating to metal sites highlight that simple metal parameter models can provide some predictive free energy estimates for metalloprotein-ligand systems, but leave room for improvement in their ease of use, modeling of coordination sites and as a result, their accuracy.

Input Pose is Key to Performance of Free Energy Perturbation: Benchmarking with Monoacylglycerol Lipase

Free energy perturbation (FEP) methodologies have become commonplace methods for modeling potency in hit-to-lead and lead optimization stages of drug discovery. The conformational states of the initial poses of compounds for FEP+ calculations are often set up by alignment to a cocrystal structure ligand, but it is not clear if this method provides the best result for all proteins or all ligands. Not only are ligand conformational states potential variables in modeling compound potency in FEP but also the selection of crystallographic water molecules for inclusion in the FEP input structures can impact FEP models. Here, we report the results of FEP calculations using FEP+ from Schrödinger and starting from maximum common substructure alignment and docked poses generated with an array of docking methodologies. As a benchmark data set, we use monoacylglycerol lipase (MAGL), an important clinical drug target in cancer malignancy, neurological diseases, and metabolic disorders, and a set of 17 MAGL inhibitors. We found a large variation among FEP+ correlations to experimental IC50 values depending on the method used to generate the input pose and that the inclusion of ligand-based information in the docking process, with some methods, increases the correlation between FEP+ free energies and IC50 values. Upon analysis of the initial poses, we found that the differences in FEP+ correlations stemmed from rotation around a tertiary amide bond as well as translation of the compound toward the more hydrophobic side of the MAGL pocket. FEP+ estimation improved across all pose modeling methods when hydrogen bond constraint information was added. However, simple maximum common substructure alignment in the presence of all crystallographic water molecules outperformed all other methods in correlation between estimated and experimental IC50 values. Taken together, these findings suggest that pose selection and crystallographic water inclusion greatly impact how well FEP+ estimated IC50 values align with experimental IC50 values and that modelers should benchmark a few different pose generation methodologies and different water inclusion strategies for their hit-to-lead and lead optimization drug discovery projects.

MM-DRPNet: A multimodal dynamic radial partitioning network for enhanced protein-ligand binding affinity prediction

Accurate prediction of drug-target binding affinity remains a fundamental challenge in contemporary drug discovery. Despite significant advances in computational methods for protein-ligand binding affinity prediction, current approaches still face substantial limitations in prediction accuracy. Moreover, the prevalent methodologies often overlook critical three-dimensional (3D) structural information, thereby constraining their practical utility in computer-aided drug design (CADD). Here we present MM-DRPNet, a multimodal deep learning framework that enhances binding affinity prediction by integrating protein-ligand structural information with interaction features and physicochemical properties. The core innovation lies in our dynamic radial partitioning (DRP) algorithm, which adaptively segments 3D space based on complex-specific interaction patterns, surpassing traditional fixed partitioning methods in capturing spatial interactions. MM-DRPNet further incorporates molecular topological features to comprehensively model both structural and spatial relationships. Extensive evaluations on benchmark datasets demonstrate that MM-DRPNet significantly outperforms state-of-the-art methods across multiple metrics, with ablation studies confirming the substantial contribution of each architectural component. Source code for MM-DRPNet is freely available for download at https://github.com/Bigrock-dd/MMDRPv1.

Superimposing Ligands with a Ligand Overlay as an Alternate Topology Model for λ-Dynamics-Based Calculations

Alchemical free energy (AFE) calculations can predict binding affinity changes as a function of structural modifications and have become powerful tools for lead optimization and drug discovery. Central to the setup and performance of AFE calculations is the manner of mapping alchemical transformations, known as the topology model. Single, dual, and hybrid topology models have been used with various AFE methods in the field. In recent works, λ-dynamics (λD) free energy calculations, specifically, have preferred the use of a hybrid multiple topology (HMT) for sampling multiple ligand perturbations. In this work, we evaluate a new topology method called ligand overlay (LO) for use with λD-based calculations, including the recently introduced λ-dynamics with a bias-updated Gibbs sampling (LaDyBUGS) approach. LO is a full multiple topology model that allows entire ligands to be sampled and restrained within a λ-dynamics framework. Relative binding free energies were computed with HMT or LO topology models with LaDyBUGS for 45 ligands across five protein benchmark systems. An overall Pearson R correlation of 0.98 and mean unsigned error of 0.32 kcal/mol were observed, suggesting that LO is a viable alternative topology model for λD-based calculations. We discuss the merits of using an HMT or LO model for future ligand studies with λD or LaDyBUGS calculations.

Combining de novo molecular design with semiempirical protein-ligand binding free energy calculation

Semi-empirical quantum chemistry methods estimate the binding free energies of protein-ligand complexes. We present an integrated approach combining the GFN2-xTB method with de novo design for the generation and evaluation of potential inhibitors of acetylcholinesterase (AChE). We employed chemical language model-based molecule generation to explore the synthetically accessible chemical space around the natural product Huperzine A, a potent AChE inhibitor. Four distinct molecular libraries were created using structure- and ligand-based molecular de novo design with SMILES and SELFIES representations, respectively. These libraries were computationally evaluated for synthesizability, novelty, and predicted biological activity. The candidate molecules were subjected to molecular docking to identify hypothetical binding poses, which were further refined using Gibbs free energy calculations. The structurally novel top-ranked molecule was chemically synthesized and biologically tested, demonstrating moderate micromolar activity against AChE. Our findings highlight the potential and certain limitations of integrating deep learning-based molecular generation with semi-empirical quantum chemistry-based activity prediction for structure-based drug design.

Binding Selectivity Analysis from Alchemical Receptor Hopping and Swapping Free Energy Calculations

We present receptor hopping and receptor swapping free energy estimation protocols based on the Alchemical Transfer Method (ATM) to model the binding selectivity of a set of ligands to two arbitrary receptors. The receptor hopping protocol, where a ligand is alchemically transferred from one receptor to another in one simulation, directly yields the ligand's binding selectivity free energy (BSFE) for the two receptors, which is the difference between the two individual binding free energies. In the receptor swapping protocol, the first ligand of a pair is transferred from one receptor to another while the second ligand is simultaneously transferred in the opposite direction. The receptor swapping free energy yields the differences in binding selectivity free energies of a set of ligands, which, when combined using a generalized DiffNet algorithm, yield the binding selectivity free energies of the ligands. We test these algorithms on host-guest systems and show that they yield results that agree with experimental data and are consistent with differences in absolute and relative binding free energies obtained by conventional methods. Preliminary applications to the selectivity analysis of molecular fragments binding to the trypsin and thrombin serine protease confirm the potential of the receptor swapping technology in structure-based drug discovery. The novel methodologies presented in this work are a first step toward streamlined and computationally efficient protocols for ligand selectivity optimization between mutants and homologous proteins.

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.

Absolute Binding Free Energies with OneOPES

The calculation of absolute binding free energies (ABFEs) for protein-ligand systems has long been a challenge. Recently, refined force fields and algorithms have improved the quality of the ABFE calculations. However, achieving the level of accuracy required to inform drug discovery efforts remains difficult. Here, we present a transferable enhanced sampling strategy to accurately calculate absolute binding free energies using OneOPES with simple geometric collective variables. We tested the strategy on two protein targets, BRD4 and Hsp90, complexed with a total of 17 chemically diverse ligands, including both molecular fragments and drug-like molecules. Our results show that OneOPES accurately predicts protein-ligand binding affinities with a mean unsigned error within 1 kcal mol-1 of experimentally determined free energies, without the need to tailor the collective variables to each system. Furthermore, our strategy effectively samples different ligand binding modes and consistently matches the experimentally determined structures regardless of the initial protein-ligand configuration. Our results suggest that the proposed OneOPES strategy can be used to inform lead optimization campaigns in drug discovery and to study protein-ligand binding and unbinding mechanisms.

Biomolecular simulations at the exascale: From drug design to organelles and beyond

The rapid advancement in computational power available for research offers to bring not only quantitative improvements, but also qualitative changes in the field of biomolecular simulation. Here, we review the state of biomolecular dynamics simulations at the threshold to exascale resources becoming available. Both developments in parallel and distributed computing will be discussed, providing a perspective on the state of the art of both. A main focus will be on obtaining binding and conformational free energies, with an outlook to macromolecular complexes and (sub)cellular assemblies.

Automated Adaptive Absolute Binding Free Energy Calculations

Alchemical absolute binding free energy (ABFE) calculations have substantial potential in drug discovery, but are often prohibitively computationally expensive. To unlock their potential, efficient automated ABFE workflows are required to reduce both computational cost and human intervention. We present a fully automated ABFE workflow based on the automated selection of λ windows, the ensemble-based detection of equilibration, and the adaptive allocation of sampling time based on inter-replicate statistics. We find that the automated selection of intermediate states with consistent overlap is rapid, robust, and simple to implement. Robust detection of equilibration is achieved with a paired t-test between the free energy estimates at initial and final portions of a an ensemble of runs. We determine reasonable default parameters for all algorithms and show that the full workflow produces equivalent results to a nonadaptive scheme over a variety of test systems, while often accelerating equilibration. Our complete workflow is implemented in the open-source package A3FE (https://github.com/michellab/a3fe).

Impact of protein conformations on binding free energy calculations in the beta-secretase 1 system

In binding free energy calculations, simulations must sample all relevant conformations of the system in order to obtain unbiased results. For instance, different ligands can bind to different metastable states of a protein, and if these protein conformational changes are not sampled in relative binding free energy calculations, the contribution of these states to binding is not accounted for and thus calculated binding free energies are inaccurate. In this work, we investigate the impact of different beta-sectretase 1 (BACE1) protein conformations obtained from x-ray crystallography on the binding of BACE1 inhibitors. We highlight how these conformational changes are not adequately sampled in typical molecular dynamics simulations. Furthermore, we show that insufficient sampling of relevant conformations induces substantial error in relative binding free energy calculations, as judged by a variation in calculated relative binding free energies up to 2 kcal/mol depending on the starting protein conformation. These results emphasize the importance of protein conformational sampling and pose this BACE1 system as a challenge case for further method development in the area of enhanced protein conformational sampling, either in combination with binding calculations or as an endpoint correction.

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.

NPEX: Never give up protein exploration with deep reinforcement learning

Elucidating unknown structures of proteins, such as metastable states, is critical in designing therapeutic agents. Protein structure exploration has been performed using advanced computational methods, especially molecular dynamics and Markov chain Monte Carlo simulations, which require untenably long calculation times and prior structural knowledge. Here, we developed an innovative method for protein structure determination called never give up protein exploration (NPEX) with deep reinforcement learning. The NPEX method leverages the soft actor-critic algorithm and the intrinsic reward system, effectively adding a bias potential without the need for prior knowledge. To demonstrate the method's effectiveness, we applied it to four models: a double well, a triple well, the alanine dipeptide, and the tryptophan cage. Compared with Markov chain Monte Carlo simulations, NPEX had markedly greater sampling efficiency. The significantly enhanced computational efficiency and lack of prior domain knowledge requirements of the NPEX method will revolutionize protein structure exploration.

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.

Robust Prediction of Relative Binding Energies for Protein-Protein Complex Mutations Using Free Energy Perturbation Calculations

Computational free energy-based methods have the potential to significantly improve throughput and decrease costs of protein design efforts. Such methods must reach a high level of reliability, accuracy, and automation to be effectively deployed in practical industrial settings in a way that impacts protein design projects. Here, we present a benchmark study for the calculation of relative changes in protein-protein binding affinity for single point mutations across a variety of systems from the literature, using free energy perturbation (FEP+) calculations. We describe a method for robust treatment of alternate protonation states for titratable amino acids, which yields improved correlation with and reduced error compared to experimental binding free energies. Following careful analysis of the largest outlier cases in our dataset, we assess limitations of the default FEP+ protocols and introduce an automated script which identifies probable outlier cases that may require additional scrutiny and calculates an empirical correction for a subset of charge-related outliers. Through a series of three additional case study systems, we discuss how Protein FEP+ can be applied to real-world protein design projects, and suggest areas of further study.

Machine-learned molecular mechanics force fields from large-scale quantum chemical data

The development of reliable and extensible molecular mechanics (MM) force fields-fast, empirical models characterizing the potential energy surface of molecular systems-is indispensable for biomolecular simulation and computer-aided drug design. Here, we introduce a generalized and extensible machine-learned MM force field, espaloma-0.3, and an end-to-end differentiable framework using graph neural networks to overcome the limitations of traditional rule-based methods. Trained in a single GPU-day to fit a large and diverse quantum chemical dataset of over 1.1 M energy and force calculations, espaloma-0.3 reproduces quantum chemical energetic properties of chemical domains highly relevant to drug discovery, including small molecules, peptides, and nucleic acids. Moreover, this force field maintains the quantum chemical energy-minimized geometries of small molecules and preserves the condensed phase properties of peptides and folded proteins, self-consistently parametrizing proteins and ligands to produce stable simulations leading to highly accurate predictions of binding free energies. This methodology demonstrates significant promise as a path forward for systematically building more accurate force fields that are easily extensible to new chemical domains of interest.

Chemical accuracy for ligand-receptor binding Gibbs energies through multi-level SQM/QM calculations

Calculating the Gibbs energies of binding of ligand-receptor systems with a thermochemical accuracy of ± 1 kcal mol-1 is a challenge to computational approaches. After exploration of the conformational space of the host, ligand and their resulting complexes upon coordination by semi-empirical GFN2 MD and meta-MD simulations, the systematic refinement through a multi-level improvement of binding modes in terms of electronic energies and solvation is able to give Gibbs energies of binding of drug molecules to CB[8] and β-CD macrocyclic receptors with such an accuracy. The accurate treatment of a small number of structures outperforms system-specific force-matching and alchemical transfer model approaches without an extensive sampling and integration.

Identification of novel LRRK2 inhibitors by structure-based virtual screening and alchemical free energy calculation

The Leucine-rich repeat kinase 2 (LRRK2) target has been identified as a promising drug target for Parkinson's disease (PD) treatment. This study focuses on optimizing the activity of LRRK2 inhibitors using alchemical relative binding free energy (RBFE) calculations. Initially, we assessed various free energy calculation methods across different LRRK2 kinase inhibitor scaffolds. The results indicate that alchemical free energy calculations are promising for prospective predictions on LRRK2 inhibitors, especially for the aminopyrimidine scaffold with an RMSE of 1.15 kcal mol-1 and Rp of 0.83. Following this, we optimized a potent LRRK2 kinase inhibitor identified from previous virtual screenings, featuring a novel scaffold. Guided by RBFE predictions using alchemical methods, this optimization led to the discovery of compound LY2023-001. This compound, with a [1,2,4]triazolo[5,6-b]indole scaffold, exhibited enhanced inhibitory activity against G2019S LRRK2 (IC50 = 12.9 nM). Molecular dynamics (MD) simulations revealed that LY2023-001 formed stable hydrogen bonds with Glu1948, and Ala1950 in the G2019S LRRK2 protein. Additionally, its phenyl substituents engage in strong electrostatic interactions with Lys1906 and van der Waals interactions with Leu1885, Phe1890, Val1893, Ile1933, Met1947, Leu1949, Leu2001, Ala2016, and Asp2017. Our findings underscore the potential of computational methods in the successful optimization of small molecules, offering important insights for the development of novel LRRK2 inhibitors.

How to Sample Dozens of Substitutions per Site with λ Dynamics

Alchemical free energy methods are useful in computer-aided drug design and computational protein design because they provide rigorous statistical mechanics-based estimates of free energy differences from molecular dynamics simulations. λ dynamics is a free energy method with the ability to characterize combinatorial chemical spaces spanning thousands of related systems within a single simulation, which gives it a distinct advantage over other alchemical free energy methods that are mostly limited to pairwise comparisons. Recently developed methods have improved the scalability of λ dynamics to perturbations at many sites; however, the size of chemical space that can be explored at each individual site has previously been limited to fewer than ten substituents. As the number of substituents increases, the volume of alchemical space corresponding to nonphysical alchemical intermediates grows exponentially relative to the size corresponding to the physical states of interest. Beyond nine substituents, λ dynamics simulations become lost in an alchemical morass of intermediate states. In this work, we introduce new biasing potentials that circumvent excessive sampling of intermediate states by favoring sampling of physical end points relative to alchemical intermediates. Additionally, we present a more scalable adaptive landscape flattening algorithm for these larger alchemical spaces. Finally, we show that this potential enables more efficient sampling in both protein and drug design test systems with up to 24 substituents per site, enabling, for the first time, simultaneous simulation of all 20 amino acids.

Rapid, Accurate, Ranking of Protein-Ligand Binding Affinities with VM2, the Second-Generation Mining Minima Method

The structure-based technologies most widely used to rank the affinities of candidate small molecule drugs for proteins range from faster but less reliable docking methods to slower but more accurate explicit solvent free energy methods. In recent years, we have advanced another technology, which is called mining minima because it "mines" out the main contributions to the chemical potentials of the free and bound molecular species by identifying and characterizing their main local energy minima. The present study provides systematic benchmarks of the accuracy and computational speed of mining minima, as implemented in the VeraChem Mining Minima Generation 2 (VM2) code, across two well-regarded protein-ligand benchmark data sets, for which there are already benchmark data for docking, free energy, and other computational methods. A core result is that VM2's accuracy approaches that of explicit solvent free energy methods at a far lower computational cost. In finer-grained analyses, we also examine the influence of various run settings, such as the treatment of crystallographic water molecules, on the accuracy, and define the costs in time and dollars of representative runs on Amazon Web Services (AWS) compute instances with various CPU and GPU combinations. We also use the benchmark data to determine the importance of VM2's correction from generalized Born to finite-difference Poisson-Boltzmann results for each energy well and find that this correction affords a remarkably consistent improvement in accuracy at a modest computational cost. The present results establish VM2 as a distinctive technology for early-stage drug discovery, which provides a strong combination of efficiency and predictivity.

A λ-Dynamics Investigation of Insulin Wakayama and Other A3 Variant Binding Affinities to the Insulin Receptor

Insulin Wakayama is a clinical insulin variant where a conserved valine at the third residue on insulin's A chain (ValA3) is replaced with a leucine (LeuA3), weakening insulin receptor (IR) binding by 140-500-fold. This severe impact on binding from a subtle modification has posed an intriguing problem for decades. Although experimental investigations of natural and unnatural A3 mutations have highlighted the sensitivity of insulin-IR binding at this site, atomistic explanations of these binding trends have remained elusive. We investigate this problem computationally using λ-dynamics free energy calculations to model structural changes in response to perturbations of the ValA3 side chain and to calculate associated relative changes in binding free energy (ΔΔGbind). The Wakayama LeuA3 mutation and seven other A3 substitutions were studied in this work. The calculated ΔΔGbind results showed high agreement compared to experimental binding potencies with a Pearson correlation of 0.88 and a mean unsigned error of 0.68 kcal/mol. Extensive structural analyses of λ-dynamics trajectories revealed that critical interactions were disrupted between insulin and the insulin receptor as a result of the A3 mutations. This investigation also quantifies the effect that adding an A3 Cδ atom or losing an A3 Cγ atom has on insulin's binding affinity to the IR. Thus, λ-dynamics was able to successfully model the effects of mutations to insulin's A3 side chain on its protein-protein interactions with the IR and shed new light on a decades-old mystery: the exquisite sensitivity of hormone-receptor binding to a subtle modification of an invariant insulin residue.

Insights and Challenges in Correcting Force Field Based Solvation Free Energies Using a Neural Network Potential

We present a comprehensive study investigating the potential gain in accuracy for calculating absolute solvation free energies (ASFE) using a neural network potential to describe the intramolecular energy of the solute. We calculated the ASFE for most compounds from the FreeSolv database using the Open Force Field (OpenFF) and compared them to earlier results obtained with the CHARMM General Force Field (CGenFF). By applying a nonequilibrium (NEQ) switching approach between the molecular mechanics (MM) description (either OpenFF or CGenFF) and the neural net potential (NNP)/MM level of theory (using ANI-2x as the NNP potential), we attempted to improve the accuracy of the calculated ASFEs. The predictive performance of the results did not change when this approach was applied to all 589 small molecules in the FreeSolv database that ANI-2x can describe. When selecting a subset of 156 molecules, focusing on compounds where the force fields performed poorly, we saw a slight improvement in the root-mean-square error (RMSE) and mean absolute error (MAE). The majority of our calculations utilized unidirectional NEQ protocols based on Jarzynski's equation. Additionally, we conducted bidirectional NEQ switching for a subset of 156 solutes. Notably, only a small fraction (10 out of 156) exhibited statistically significant discrepancies between unidirectional and bidirectional NEQ switching free energy estimates.

MDFit: automated molecular simulations workflow enables high throughput assessment of ligands-protein dynamics

Molecular dynamics (MD) simulation is a powerful tool for characterizing ligand-protein conformational dynamics and offers significant advantages over docking and other rigid structure-based computational methods. However, setting up, running, and analyzing MD simulations continues to be a multi-step process making it cumbersome to assess a library of ligands in a protein binding pocket using MD. We present an automated workflow that streamlines setting up, running, and analyzing Desmond MD simulations for protein-ligand complexes using machine learning (ML) models. The workflow takes a library of pre-docked ligands and a prepared protein structure as input, sets up and runs MD with each protein-ligand complex, and generates simulation fingerprints for each ligand. Simulation fingerprints (SimFP) capture protein-ligand compatibility, including stability of different ligand-pocket interactions and other useful metrics that enable easy rank-ordering of the ligand library for pocket optimization. SimFPs from a ligand library are used to build & deploy ML models that predict binding assay outcomes and automatically infer important interactions. Unlike relative free-energy methods that are constrained to assess ligands with high chemical similarity, ML models based on SimFPs can accommodate diverse ligand sets. We present two case studies on how SimFP helps delineate structure-activity relationship (SAR) trends and explain potency differences across matched-molecular pairs of (1) cyclic peptides targeting PD-L1 and (2) small molecule inhibitors targeting CDK9.

Current State of Open Source Force Fields in Protein-Ligand Binding Affinity Predictions

In drug discovery, the in silico prediction of binding affinity is one of the major means to prioritize compounds for synthesis. Alchemical relative binding free energy (RBFE) calculations based on molecular dynamics (MD) simulations are nowadays a popular approach for the accurate affinity ranking of compounds. MD simulations rely on empirical force field parameters, which strongly influence the accuracy of the predicted affinities. Here, we evaluate the ability of six different small-molecule force fields to predict experimental protein-ligand binding affinities in RBFE calculations on a set of 598 ligands and 22 protein targets. The public force fields OpenFF Parsley and Sage, GAFF, and CGenFF show comparable accuracy, while OPLS3e is significantly more accurate. However, a consensus approach using Sage, GAFF, and CGenFF leads to accuracy comparable to OPLS3e. While Parsley and Sage are performing comparably based on aggregated statistics across the whole dataset, there are differences in terms of outliers. Analysis of the force field reveals that improved parameters lead to significant improvement in the accuracy of affinity predictions on subsets of the dataset involving those parameters. Lower accuracy can not only be attributed to the force field parameters but is also dependent on input preparation and sampling convergence of the calculations. Especially large perturbations and nonconverged simulations lead to less accurate predictions. The input structures, Gromacs force field files, as well as the analysis Python notebooks are available on GitHub.

Complexation Mechanisms of Aqueous Amylose: Molecular Dynamics Study Using 3-Pentadecylphenol

Molecular dynamics (MD) simulations of linear amylose fragments containing 10 to 40 glucose units were used to study the complexation of the prototypical compound, 3-pentadecylphenol (PDP)─a natural product with surfactant-like properties─in aqueous solution. The amylose-PDP binding leverages mainly hydrophobic interactions together with excluded volume effects. It was found that while the most stable complexes contained PDP inside the helical structure of the amylose in the expected guest-host (inclusion) complexation manner, at higher temperatures, the commonly observed PDP-amylose complexes often involved more nonspecific interactions than inclusion complexation. In the case where a stoichiometric excess of PDP was added to the simulation box, self-aggregation of the small molecule precluded its ability to enter the internal helical part of the oligosaccharide, and as a result, inclusion complexation became ineffective. MD simulation trajectories were analyzed preliminarily using cluster analysis (CA), followed by more rigorous solvent accessible surface area (SASA) determination over the temperature range spanning from 277 to 433 K. It was found that using the SASA of PDP corrected for its intrinsic conformational changes, together with a generic hidden Markov model (HMM), an adequate quantification of the different types of PDP-amylose aggregates was obtained to allow further analysis. The enthalpy change associated with the guest-host binding equilibrium constant (Kgh) in aqueous solution was estimated to be -75 kJ/mol, which is about twice as high as one might expect based on experimentally measured values of similar complexes in the solid state where the (unsolvated) helical structure of amylose remains rigid. On the other hand, the nonspecific binding (Kns) enthalpy change associated with PDP-amylose interactions in the same solution environment was found to be about half of the inclusion complexation value.

DrugGym: A testbed for the economics of autonomous drug discovery

Drug discovery is stochastic. The effectiveness of candidate compounds in satisfying design objectives is unknown ahead of time, and the tools used for prioritization-predictive models and assays-are inaccurate and noisy. In a typical discovery campaign, thousands of compounds may be synthesized and tested before design objectives are achieved, with many others ideated but deprioritized. These challenges are well-documented, but assessing potential remedies has been difficult. We introduce DrugGym, a framework for modeling the stochastic process of drug discovery. Emulating biochemical assays with realistic surrogate models, we simulate the progression from weak hits to sub-micromolar leads with viable ADME. We use this testbed to examine how different ideation, scoring, and decision-making strategies impact statistical measures of utility, such as the probability of program success within predefined budgets and the expected costs to achieve target candidate profile (TCP) goals. We also assess the influence of affinity model inaccuracy, chemical creativity, batch size, and multi-step reasoning. Our findings suggest that reducing affinity model inaccuracy from 2 to 0.5 pIC50 units improves budget-constrained success rates tenfold. DrugGym represents a realistic testbed for machine learning methods applied to the hit-to-lead phase. Source code is available at www.drug-gym.org.

Multistate Method to Efficiently Account for Tautomerism and Protonation in Alchemical Free-Energy Calculations

The majority of drug-like molecules contain at least one ionizable group, and many common drug scaffolds are subject to tautomeric equilibria. Thus, these compounds are found in a mixture of protonation and/or tautomeric states at physiological pH. Intrinsically, standard classical molecular dynamics (MD) simulations cannot describe such equilibria between states, which negatively impacts the prediction of key molecular properties in silico. Following the formalism described by de Oliveira and co-workers (J. Chem. Theory Comput. 2019, 15, 424-435) to consider the influence of all states on the binding process based on alchemical free-energy calculations, we demonstrate in this work that the multistate method replica-exchange enveloping distribution sampling (RE-EDS) is well suited to describe molecules with multiple protonation and/or tautomeric states in a single simulation. We apply our methodology to a series of eight inhibitors of factor Xa with two protonation states and a series of eight inhibitors of glycogen synthase kinase 3β (GSK3β) with two tautomeric states. In particular, we show that given a sufficient phase-space overlap between the states, RE-EDS is computationally more efficient than standard pairwise free-energy methods.

Robust prediction of relative binding energies for protein-protein complex mutations using free energy perturbation calculations

Computational free energy-based methods have the potential to significantly improve throughput and decrease costs of protein design efforts. Such methods must reach a high level of reliability, accuracy, and automation to be effectively deployed in practical industrial settings in a way that impacts protein design projects. Here, we present a benchmark study for the calculation of relative changes in protein-protein binding affinity for single point mutations across a variety of systems from the literature, using free energy perturbation (FEP+) calculations. We describe a method for robust treatment of alternate protonation states for titratable amino acids, which yields improved correlation with and reduced error compared to experimental binding free energies. Following careful analysis of the largest outlier cases in our dataset, we assess limitations of the default FEP+ protocols and introduce an automated script which identifies probable outlier cases that may require additional scrutiny and calculates an empirical correction for a subset of charge-related outliers. Through a series of three additional case study systems, we discuss how protein FEP+ can be applied to real-world protein design projects, and suggest areas of further study.

Prospective de novo drug design with deep interactome learning

De novo drug design aims to generate molecules from scratch that possess specific chemical and pharmacological properties. We present a computational approach utilizing interactome-based deep learning for ligand- and structure-based generation of drug-like molecules. This method capitalizes on the unique strengths of both graph neural networks and chemical language models, offering an alternative to the need for application-specific reinforcement, transfer, or few-shot learning. It enables the "zero-shot" construction of compound libraries tailored to possess specific bioactivity, synthesizability, and structural novelty. In order to proactively evaluate the deep interactome learning framework for protein structure-based drug design, potential new ligands targeting the binding site of the human peroxisome proliferator-activated receptor (PPAR) subtype gamma are generated. The top-ranking designs are chemically synthesized and computationally, biophysically, and biochemically characterized. Potent PPAR partial agonists are identified, demonstrating favorable activity and the desired selectivity profiles for both nuclear receptors and off-target interactions. Crystal structure determination of the ligand-receptor complex confirms the anticipated binding mode. This successful outcome positively advocates interactome-based de novo design for application in bioorganic and medicinal chemistry, enabling the creation of innovative bioactive molecules.

AI is a viable alternative to high throughput screening: a 318-target study

High throughput screening (HTS) is routinely used to identify bioactive small molecules. This requires physical compounds, which limits coverage of accessible chemical space. Computational approaches combined with vast on-demand chemical libraries can access far greater chemical space, provided that the predictive accuracy is sufficient to identify useful molecules. Through the largest and most diverse virtual HTS campaign reported to date, comprising 318 individual projects, we demonstrate that our AtomNet® convolutional neural network successfully finds novel hits across every major therapeutic area and protein class. We address historical limitations of computational screening by demonstrating success for target proteins without known binders, high-quality X-ray crystal structures, or manual cherry-picking of compounds. We show that the molecules selected by the AtomNet® model are novel drug-like scaffolds rather than minor modifications to known bioactive compounds. Our empirical results suggest that computational methods can substantially replace HTS as the first step of small-molecule drug discovery.

Enhancing Generalizability in Protein-Ligand Binding Affinity Prediction with Multimodal Contrastive Learning

Improving the generalization ability of scoring functions remains a major challenge in protein-ligand binding affinity prediction. Many machine learning methods are limited by their reliance on single-modal representations, hindering a comprehensive understanding of protein-ligand interactions. We introduce a graph-neural-network-based scoring function that utilizes a triplet contrastive learning loss to improve protein-ligand representations. In this model, three-dimensional complex representations and the fusion of two-dimensional ligand and coarse-grained pocket representations converge while distancing from decoy representations in latent space. After rigorous validation on multiple external data sets, our model exhibits commendable generalization capabilities compared to those of other deep learning-based scoring functions, marking it as a promising tool in the realm of drug discovery. In the future, our training framework can be extended to other biophysical- and biochemical-related problems such as protein-protein interaction and protein mutation prediction.