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

Discovery of Novel Multiangiogenic Agents Targeting VEGFR2, EphB4, FGFR-1, and TIE-2: Receptor-Based Pharmacophore Modeling, Virtual Screening, and Molecular Modeling Studies

The angiogenesis phenomenon is crucial for the formation of new blood vessels in cancer cells. The cancerous cells' progress hampers other healthy cells. The main objective of this study is to explore and decipher multimodal natural compounds against VEGFR2, EphB4, FGFR-1, and TIE-2 drug targets to arrest angiogenesis and progression. The receptor-based pharmacophore modeling of VEGFR2, EphB4, FGFR-1, and TIE-2 was developed and validated through enrichment parameters. Further, the validated hypothesis allowed for screening druglike natural product databases such as SuperNatural 3.0, COCONUT, and LOTUS. The common pharmacophoric featured natural compounds were assessed for binding affinities using absolute end-point methods. Finally, density functional theory has been studied to understand the chemical reactivity and stability of the protein complexes. Among all of the screened natural compounds, 17 natural compounds were found to align accurately against validated pharmacophore models having higher fitness scores and align scores. Taking reference drugs sorafenib (VEGFR2), NVP-BHG712 (EphB4), pemiganitib (FGFR-1), and DP1919 (TIE-2), three promising natural compounds CNP0003920, CNP0243075, and CNP0211397 were concluded based on their end-point binding energies, binding interactions, molecular dynamics, and optimal pharmacokinetic and toxicity profiles. The density functional theory (DFT) results suggested that the identified compounds bound with protein complexes are stable. Our findings can represent a promising starting point for developing multimodal analogues VEGFR2, EphB4, FGFR-1, and TIE-2 proteins.

Divide-and-Conquer ABFE: Improving Free Energy Calculations by Enhancing Water Sampling

Free energy perturbation (FEP) is a promising method for accurately predicting molecular interactions, widely applied in fields such as drug design, materials science, and catalysis. However, FEP calculations, particularly those in aqueous environments, often suffer from convergence issues due to insufficient sampling of water molecules. These challenges are particularly critical in solvation-related free energy calculations, such as small molecule-protein binding, interface interactions, and molecular adsorption on surfaces. To address these limitations, we present the divide-and-conquer absolute binding free energy (DC-ABFE) method. By partitioning the ligand or molecule into atomic groups and sequentially decoupling their van der Waals interactions, DC-ABFE improves water re-entry sampling, enhances phase-space overlap, and significantly enhances the convergence of free energy calculations. Our benchmark demonstrates that DC-ABFE achieves more reproducible and reliable binding free energy predictions compared to traditional FEP methods. DC-ABFE is applicable to a range of free energy calculations involving solvation effects. Additionally, this method establishes a stronger foundation for precise drug screening, offering a more robust tool for predicting binding affinities in drug discovery.

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.

Robust Automated Truncation Point Selection for Molecular Simulations

Quantities calculated from molecular simulations are often subject to an initial bias due to unrepresentative starting configurations. Initial data are usually discarded to reduce bias. Chodera's method for automated truncation point selection [J. Chem. Theory Comput. 2016, 12, 4, 1799-1805] is popular but has not been thoroughly assessed. We reformulate White's marginal standard error rule to provide a spectrum of truncation point selection heuristics that differ in their treatment of autocorrelation. These include a method effectively equivalent to Chodera's. We test these methods on ensembles of synthetic time series modeled on free energy change estimates from long absolute binding free energy calculations. Methods that more thoroughly account for autocorrelation often show late and variable truncation times, while methods that less thoroughly account for autocorrelation often show early truncation, relative to the optimal truncation point. This increases variance and bias, respectively. We recommend a method that achieves robust performance across our test sets by balancing these two extremes. None of the methods reliably detected insufficient sampling. All heuristics tested are implemented in the open-source Python package RED (github.com/fjclark/red).

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.

In silico studies for improving target selectivity of anti-malarial dual falcipain inhibitors vis-à-vis human cathepsins

Dual falcipain-2 (FP-2) and falcipain-3 (FP-3) inhibitors, NM12 and NM15, displayed micromolar inhibitions but they exhibit similar binding affinities for the human cathepsins, thus indicating potential toxicity. The current study aims to develop a model to enhance the selectivity of the falcipain inhibitors vis-à-vis human cathepsins using previously identified dual falcipain 2 and 3 inhibitors, NM12 and NM15. To improve the selectivity of NM12 and NM15, analogs with weaker interactions with the conserved residues in the FPs and hCatK were designed while enhancing the unique interactions for the FPs. In silico analysis was carried out in the S2 subsite of both plasmodium and human proteases which is considered the preferred selective site due to the presence of less conserved residues. The Fasta sequence alignment and active/conserved binding site superimposition show that FPs contain acidic polar residues (Asp234 for FP2 and Glu243 for FP3) while hCatK has a neutral hydrophobic residue (Leu209) at the S2 subsite. Therefore, analogs of NM12 and NM15 were designed to enhance affinity and selectivity by improving interactions with these acidic residues while avoiding interactions with hydrophobic residues in hCatK. Newly designed analogs (NM12H and NM15G) show better selectivity as well as binding affinity towards FPs (ΔG of NM12H: -74.49 kcal/mol for FP2, -70.97 kcal/mol for FP3; ΔG of NM15G: -70.09 kcal/mol for FP2, -74.52 kcal/mol for FP3) as compared to NM12 and NM15. Thus, the selectivity and binding affinity against dual falcipains vis-à-vis human cathepsin were improved using molecular dynamic simulations.

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.

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

Lambda-ABF: Simplified, Portable, Accurate, and Cost-Effective Alchemical Free-Energy Computation

We introduce the lambda-Adaptive Biasing Force (lambda-ABF) method for the computation of alchemical free-energy differences. We propose a software implementation and showcase it on biomolecular systems. The method arises from coupling multiple-walker adaptive biasing force with λ-dynamics. The sampling of the alchemical variable is continuous and converges toward a uniform distribution, making manual optimization of the λ schedule unnecessary. Contrary to most other approaches, alchemical free-energy estimates are obtained immediately without any postprocessing. Free diffusion of λ improves orthogonal relaxation compared to fixed-λ thermodynamic integration or free-energy perturbation. Furthermore, multiple walkers provide generic orthogonal space coverage with minimal user input and negligible computational overhead. We show that our high-performance implementations coupling the Colvars library with NAMD and Tinker-HP can address real-world cases including ligand-receptor binding with both fixed-charge and polarizable models, with a demonstrably richer sampling than fixed-λ methods. The implementation is fully open-source, publicly available, and readily usable by practitioners of current alchemical methods. Thanks to the portable Colvars library, lambda-ABF presents a unified user interface regardless of the back-end (NAMD, Tinker-HP, or any software to be interfaced in the future), sparing users the effort of learning multiple interfaces. Finally, the Colvars Dashboard extension of the visual molecular dynamics (VMD) software provides an interactive monitoring and diagnostic tool for lambda-ABF simulations.

Sire: An interoperability engine for prototyping algorithms and exchanging information between molecular simulation programs

Sire is a Python/C++ library that is used both to prototype new algorithms and as an interoperability engine for exchanging information between molecular simulation programs. It provides a collection of file parsers and information converters that together make it easier to combine and leverage the functionality of many other programs and libraries. This empowers researchers to use sire to write a single script that can, for example, load a molecule from a PDBx/mmCIF file via Gemmi, perform SMARTS searches via RDKit, parameterize molecules using BioSimSpace, run GPU-accelerated molecular dynamics via OpenMM, and then display the resulting dynamics trajectory in a NGLView Jupyter notebook 3D molecular viewer. This functionality is built on by BioSimSpace, which uses sire's molecular information engine to interconvert with programs such as GROMACS, NAMD, Amber, and AmberTools for automated molecular parameterization and the running of molecular dynamics, metadynamics, and alchemical free energy workflows. Sire comes complete with a powerful molecular information search engine, plus trajectory loading and editing, analysis, and energy evaluation engines. This, when combined with an in-built computer algebra system, gives substantial flexibility to researchers to load, search for, edit, and combine molecular information from multiple sources and use that to drive novel algorithms by combining functionality from other programs. Sire is open source (GPL3) and is available via conda and at a free Jupyter notebook server at https://try.openbiosim.org. Sire is supported by the not-for-profit OpenBioSim community interest company.

On Analytical Corrections for Restraints in Absolute Binding Free Energy Calculations

Double decoupling absolute binding free energy simulations require an intermediate state at which the ligand is held solely by restraints in a position and orientation resembling the bound state. One possible choice consists of one distance, two angle, and three dihedral angle restraints. Here, I demonstrate that in practically all cases the analytical correction derived under the rigid rotator harmonic oscillator approximation is sufficient to account for the free energy of the restraints.

Nonbonded Force Field Parameters from MBIS Partitioning of the Molecular Electron Density Improve Binding Affinity Predictions of the T4-Lysozyme Double Mutant

The use of computer simulation for binding affinity prediction is growing in drug discovery. However, its wider use is constrained by the accuracy of the free energy calculations. The key sources of error are the force fields used to depict molecular interactions and insufficient sampling of the configurational space. To improve the quality of the force field, we developed a Python-based computational workflow. The workflow described here uses the minimal basis iterative stockholder (MBIS) method to determine atomic charges and Lennard-Jones parameters from the polarized molecular density. This is done by performing electronic structure calculations on various configurations of the ligand when it is both bound and unbound. In addition, we validated a simulation procedure that accounts for the protein and ligand degrees of freedom to precisely calculate binding free energies. This was achieved by comparing the self-adjusted mixture sampling and nonequilibrium thermodynamic integration methods using various protein and ligand conformations. The accuracy of predicting binding affinity is improved by using MBIS-derived force field parameters and a validated simulation procedure. This improvement surpasses the chemical precision for the eight aromatic ligands, reaching a root-mean-square error of 0.7 kcal/mol.

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.

Tunnel-like region observed as a potential allosteric site in Staphylococcus aureus Glyceraldehyde-3-phosphate dehydrogenase

Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) catalyzing the sixth step of glycolysis has been investigated for allosteric features that might be used as potential target for specific inhibition of Staphylococcus aureus (S.aureus). X-ray structure of bacterial enzyme for which a tunnel-like opening passing through the center previously proposed as an allosteric site has been subjected to six independent 500 ns long Molecular Dynamics simulations. Harmonic bond restraints were employed at key residues to underline the allosteric feature of this region. A noticeable reduction was observed in the mobility of NAD+ binding domains when restrictions were applied. Also, a substantial decrease in cross-correlations between distant Cα fluctuations was detected throughout the structure. Mutual information (MI) analysis revealed a similar decrease in the degree of correspondence in positional fluctuations in all directions everywhere in the receptor. MI between backbone and side-chain torsional variations changed its distribution profile and decreased considerably around the catalytic sites when restraints were employed. Principal component analysis clearly showed that the restrained state sampled a narrower range of conformations than apo state, especially in the first principal mode due to restriction in the conformational flexibility of NAD+ binding domain. Clustering the trajectory based on catalytic site residues displayed a smaller repertoire of conformations for restrained state compared to apo. Representative snapshots subjected to k-shortest pathway analysis revealed the impact of bond restraints on the allosteric communication which displayed distinct optimal and suboptimal pathways for two states, where observed frequencies of critical residues Gln51 and Val283 at the proposed site changed considerably.

Accelerating and Automating the Free Energy Perturbation Absolute Binding Free Energy Calculation with the RED-E Function

The accurate prediction of the binding affinities between small molecules and biological macromolecules plays a fundamental role in structure-based drug design, which is still challenging. The free energy perturbation-based absolute binding free energy (FEP-ABFE) approach has shown potential in its reliability. To correctly calculate the energy related to the ligand being restrained by the receptor, additional restraints between the ligand and the receptor are needed. However, determining the restraint parameters for individual ligands empirically is too trivial to be automated, and usually gives rise to numerical instabilities, which set back the applications of FEP-ABFE. To address these issues, we derived the analytical expression for the probability distribution of energy differences, P(ΔU), during the process of restraint addition, which is called the RED-E (restraint energy distribution at equilibrium position) function. Simulations indicated that the RED-E function can accurately describe P(ΔU) when restraints are added at the equilibrium position. Based on the RED-E function, an automatic restraint selection method was proposed to select the best restraint. With this method, there is a high phase-space overlap between the free and restrained states, such that using a 2-λ perturbation can accurately calculate the free energy of the restraint addition, which is a nearly 6 times acceleration compared with current widely used 12-λ perturbation method. The RED-E function gives insight into the non-Gaussian behavior of the sampled P(ΔU) in certain FEP processes in an analytical way. The highly automated and accelerated restraint selection also makes it possible for the large-scale application of FEP-ABFE in real drug discovery practices.

Expanded ensemble predictions of absolute binding free energies in the SAMPL9 host-guest challenge

As part of the SAMPL9 community-wide blind host-guest challenge, we implemented an expanded ensemble workflow to predict absolute binding free energies for 13 small molecules against pillar[6]arene. Notable features of our protocol include consideration of a variety of protonation and enantiomeric states for both host and guests, optimization of alchemical intermediates, and analysis of free energy estimates and their uncertainty using large numbers of simulation replicates performed using distributed computing. Our predictions of absolute binding free energies resulted in a mean absolute error of 2.29 kcal mol-1 and an R2 of 0.54. Overall, results show that expanded ensemble calculations using all-atom molecular dynamics simulations are a valuable and efficient computational tool in predicting absolute binding free energies.

Alchemical Free Energy Calculations on Membrane-Associated Proteins

Membrane proteins have diverse functions within cells and are well-established drug targets. The advances in membrane protein structural biology have revealed drug and lipid binding sites on membrane proteins, while computational methods such as molecular simulations can resolve the thermodynamic basis of these interactions. Particularly, alchemical free energy calculations have shown promise in the calculation of reliable and reproducible binding free energies of protein-ligand and protein-lipid complexes in membrane-associated systems. In this review, we present an overview of representative alchemical free energy studies on G-protein-coupled receptors, ion channels, transporters as well as protein-lipid interactions, with emphasis on best practices and critical aspects of running these simulations. Additionally, we analyze challenges and successes when running alchemical free energy calculations on membrane-associated proteins. Finally, we highlight the value of alchemical free energy calculations calculations in drug discovery and their applicability in the pharmaceutical industry.