Viral proteases are central targets in antiviral drug discovery and development because they play essential roles in viral replication and maturation. Although protease inhibitors have achieved major clinical success, traditional design strategies face challenges, including resistance development, poor oral exposure of early peptidomimetics, and off-target toxicity of highly reactive covalent warheads. Classical approaches, such as peptidomimetics, macrocyclization, and covalent warhead engineering, are discussed alongside contemporary strategies, including allosteric modulation and targeted protease degradation via proteolysis-targeting chimeras (PROTAC) technology. Particular emphasis is placed on how these strategies address key obstacles, such as resistance evolution, selectivity, metabolic stability, and oral bioavailability. Several quantitative case studies have also demonstrated the growing significance of computational tools in contemporary antiviral discovery. For SARS-CoV-2 main protease (Mpro), these workflows were enabled by the rapid availability of high-resolution experimental crystal structures of the target protein. The evolution of a weak fragment (Kd ≈ 1.7 mM; ΔG ≈ -3.6 kcal/mol) into a covalent inhibitor (QUB-00006-Int-07) with enzymatic inhibition (IC50 ≈ 830 nM) was successfully guided by molecular dynamics (MD) simulations and absolute binding free energy calculations. This was subsequently confirmed experimentally using NMR, ESI-MS, and FRET assays. Furthermore, out of 25 computationally prioritized candidates with Ki values less than 4 μM, 15 active Mpro inhibitors were identified using accelerated free-energy perturbation-based repurposing campaigns. Long-range allosteric pathways connecting the catalytic site to resistance-associated regions and experimentally verified allosteric pockets have also been discovered using dynamic nonequilibrium MD. Together, these integrated in silico approaches enable the early prioritization of high-affinity ligands, mechanistic understanding of resistance, and significant reduction of late-stage attrition in antiviral drug discovery. Through detailed case studies on SARS-CoV-2 main protease (Mpro), Zika virus NS2B-NS3 protease, and Dengue virus NS2B-NS3 protease, the review illustrates how medicinal chemistry principles translate molecular insights into clinically relevant antivirals. Finally, a forward-looking development roadmap is proposed that integrates potency, selectivity, pharmacokinetics, manufacturability, and resistance management toward the goal of broad-spectrum, durable, and adaptable protease-targeted therapeutics development.

Absolute binding free energy (ΔG) calculations can rank structurally diverse compounds, which could be useful for early-stage drug discovery. Unfortunately, for flexible systems, it can be challenging to sample the receptor conformations necessary to obtain converged ΔG calculations. Here, we address this challenge by leveraging extensive molecular dynamics simulations of apo SARS-CoV-2 main protease (MPro) that were conducted on the Folding@Home distributed computing system. A Markov state model (MSM) was built to compute the equilibrium probability of each snapshot. Representative snapshots were selected from clusters defined based on occupancy fingerprints of the catalytic site. The binding potential of mean force (BPMF), the binding free energy between a ligand and rigid receptor configuration, was computed between the representative snapshots and 130 drug leads from the COVID Moonshot, an open-source drug discovery project. ΔGs were computed using an exponential average of BPMFs based on implicit ligand theory (ILT). ΔG calculations recapitulated experimental values with a Pearson R of 0.55 and a mean-adjusted root-mean-square error of 1.6 kcal/mol. Accuracy and computational costs were found to be intermediate between docking and previous free energy calculations with a fully flexible receptor. Moreover, in 88% of systems, the calculated ΔG of the native binding pose (RMSD from crystallographic <3 Å) was within 1 kT of the top-ranked pose.

Docking-based virtual screening (VS) is essential for hit finding in the initial stage of drug or probe discovery. However, it remains prone to high false-positive rates, often resulting in unsuccessful screening campaigns. MD-based alchemical free-energy methods offer a promising solution to improve VS hit rates but are highly resource-intensive. Real-world and benchmark studies incorporating alchemical absolute binding free energy (ABFE) calculations could help optimize their use in VS pipelines. Here, we present a large-scale benchmark to evaluate the comparative value of ABFE calculations in VS workflows. Two data sets were used: a curated set of 632 ligand-protein complexes from the PDBbind database to assess ABFE quantitative accuracy and a set of 315 binders and decoys from the Database of Useful Decoys (DUD-E) to evaluate predictive power in a VS context. Alongside alchemical ABFE, we benchmarked computationally affordable end-state physics-based methods and five machine-learning (ML) models. The study ranked BFE predictors consistently with their computational cost, with alchemical ABFE performing well across both benchmarks. End-state methods scored well in recognizing actives from decoys in the DUD-E data set but showed little correlation with experimental values in PDBbind. Most ML models performed well on PDBbind, likely due to training overlap, but failed on DUD-E, except for GNINA and Boltz-2, which demonstrated a degree of generalization comparable to end-state physics-based methods. Overall, a staged approach involving Boltz-2 as a primary filter followed by alchemical ABFE is likely to robustly and cost-efficiently enrich docking-based VS hit lists with true actives.

TRPV5 is a calcium-selective epithelial ion channel that is essential for renal Ca2+ reabsorption and systemic calcium homeostasis. Despite recent structural advances, the energetic basis and coordinated contributions of the three selectivity-filter residues (D542, T539, and N572/I575) to Ca2+ selectivity remain unresolved. Here, we combine absolute and relative binding free-energy calculations (ABFE/RBFE), pore radius, and electrostatic analyses, and adaptive steered molecular dynamics (ASMD) to quantify how Ca2+ and Na+ interact with and traverse these key sites in both wild-type and mutant channels. Free-energy analyses show that all three sites favor Ca2+ over Na+, with the outer D542 site making the dominant contribution to Ca2+ selectivity, whereas T539 and N572/I575 provide weaker secondary contributions. ASMD-derived free-energy profiles further reveal a pronounced energetic preference for the Ca2+-driven displacement of Na+. Mutation of any one of the three selectivity-filter residues induced perturbations of pore radius and electrostatic landscapes attenuated the energetic preference, revealing the cooperative contribution of these three sites to Ca2+ selectivity in wild-type TRPV5. Based on the above results, we provide a hierarchical mechanism involving high-affinity Ca2+ capture at D542 and subsequent downstream modulation of the free-energy landscape by T539 and N572, forming an integrated structure-energy framework for selective Ca2+ permeation, offering testable hypotheses for disease variants and targeted modulation of TRPV5 in calcium-handling disorders.

We report a large thermodynamic effect of solvent buffer on the standard binding free energy for a large hydrophobic ligand and show that a realistic comparison with the experimental binding affinity requires correctly accounting for the solvent reference state and ligand reorganization. We focus on lenacapavir (LEN; MW ≈ 1 kDa), an HIV-1 capsid inhibitor with very low aqueous solubility. Using several absolute binding free energy (ABFE) protocols including double-decoupling method (DDM), potential-of-mean-force (PMF) approaches, and the Alchemical Transfer Method (ATM), we obtained standard binding free energy values in neat water of ΔGbind0 = -26.4 to -30.0 kcal/mol for LEN binding to the HIV-1 capsid CA dimer, much stronger than the SPR-derived affinity measured in 5% DMSO buffer (≈-13.4 kcal/mol). We analyze the discrepancy and identify two dominant contributors to the calculated overbinding. (i) Solvent reference state: A thermodynamic cycle analysis and solvation free energy calculations reveal that the 5% DMSO buffer stabilizes the free ligand by ∼-4 kcal/mol relative to neat water. Additionally, we propose that the enrichment of the hydrophobic cosolvent in the apo binding pocket makes it more energetically costly to displace DMSO upon ligand binding. (ii) Ligand reorganization: incomplete treatment of LEN's internal conformational reorganization in the ABFE protocols leads to overbinding; a DDM variant with explicit ligand reorganization reduces the overestimate by ∼6 kcal/mol. Together, considerations of these effects significantly reduce the discrepancy between the ABFE calculations and experiments. Our results suggest that for large, hydrophobic ligands, quantitative agreement between ABFEs and experiments requires (a) reporting ΔGbind0 in the appropriate assay buffer (not simply in water) and (b) explicit treatment of ligand reorganization.

In preclinical models, Phosphodiesterase-10A (PDE10A) inhibition has shown efficacy for Parkinson's and Huntington's diseases and is potentially a therapeutic approach for schizophrenia. In this study, computational approaches were employed to identify selective PDE10A inhibitors from a series of olefinated benzosuberene analogues. The molecular interactions and docking scores of the selected molecules were compared with three different co-crystal inhibitors of PDE10A. Molecular dynamics (MD) simulations confirmed the stability of the screened hit molecules with PDE10A, and AMKPD-52 was found to uniquely interact with the selectivity pocket residue Tyr-683. Further refinement using steered MD simulations and umbrella sampling simulations with a total simulation time of ~72 μs reinforced AMKPD-52 as the potential lead molecule. Absolute binding free energy calculations performed using multistate Bennett acceptance ratio method consistently ranked AMKPD-52 as the potential lead molecule for PDE10A inhibition. Experimental evaluation revealed an IC₅₀ value of 11.52 μM for AMKPD-52, confirming its inhibitory activity against PDE10A. Although the binding affinity and potency were modest compared to the reference inhibitor TAK-063, the natural origin and straightforward synthetic tractability of AMKPD-52 make it an attractive candidate for further medicinal chemistry optimization. These findings warrant in vivo validation to assess the therapeutic potential of AMKPD-52 in PDE10A-targeted interventions.

We present BindFlow, a Python-based software for automated absolute binding free energy (ABFE) calculations at the free energy perturbation (FEP) or at the molecular mechanics Poisson-Boltzmann/generalized Born surface area [MM(PB/GB)SA] level of theory. BindFlow is free, open-source, user-friendly, and easily customizable, runs on workstations or distributed computing platforms, and provides extensive documentation and tutorials. BindFlow uses GROMACS as a molecular dynamics engine and provides built-in support for the small-molecule force fields GAFF, OpenFF, and Espaloma, as well as support for user-provided custom force fields. We test BindFlow by computing affinities for 139 receptor-ligand pairs, involving eight different targets, including six soluble proteins, one membrane protein, and one nonprotein host-guest system. We find that the agreement of BindFlow predictions with experiments is overall similar to gold standards in the field. Interestingly, we find that MM(PB/GB)SA achieves correlations that, for some systems and force fields, approach those obtained with FEP while requiring only a fraction of the computational cost. This study establishes BindFlow as a validated and accessible tool for ABFE calculations.

Absolute (ABFE) and relative (RBFE) binding free energy calculations with all-atom molecular dynamics (MD) can significantly reduce costs in the early stages of drug discovery. We introduce a new implementation of the Binding Affinity Tool (BAT.py) software, which adds RBFE calculations using separated topologies (SepTop) to the already established ABFE fully automated workflow. SepTop combines the advantages of ABFE and RBFE, being applicable to ligands that have very little or no similarity, while at the same time avoiding common challenges of ABFE calculations, such as occluded binding sites and problematic conformational changes of the receptor upon ligand binding. Three different thermodynamic paths for the relative calculations were implemented into the BAT software using the AMBER and OpenMM simulation engines, and here we test them on the BRD4(2) benchmark system. We discuss their correlation with ABFE, standard RBFE, and experimental results, and also their associated computational cost.

Targeting RNA with small molecules represents a promising yet relatively unexplored avenue for the design of new drugs. Nevertheless, challenges arise from the lack of computational models and techniques able to accurately model RNA systems, and predict their binding affinities to small molecules. Here, we tackle these difficulties by developing a tailored state-of-the-art approach for absolute binding free energy calculations of RNA-binding small molecules. For this, we combine the advanced AMOEBA polarizable force field to the newly developed lambda-Adaptive Biasing Force scheme associated to refined restraints allowing for efficient sampling. To capture the free energy barrier associated to challenging RNA conformational changes, we combine machine learning-based collective variables with enhanced sampling simulations. Applying this computational protocol to a complex Riboswitch-like RNA target demonstrates quantitative predictions. These results pave the way for the routine application of free energy simulations in RNA-targeted drug discovery, thus providing a significant reduction in their failure rate.

Receptor tyrosine kinase-like orphan receptor 1 (ROR1) is a cancer-associated pseudokinase with low expression in normal adult tissues but elevated levels in various malignancies, making it a promising therapeutic target. Among ~ 4 million compounds, CHEMBL3926946 emerged as the most promising candidate, demonstrating a persistent binding pose and a well-defined free energy basin. Well-tempered metadynamics (wt-MetaD) revealed a deep minimum of 26.00 ± 2.44kcal/mol, indicating a highly stable interaction. CHEMBL3926946 exhibited a favourable Absolute Binding Free Energy Perturbation (ABFEP) of - 16.52 ± 0.37 kcal/mol, significantly outperforming the inhibitor Ponatinib (- 8.67 ± 0.94 kcal/mol), supported by persistent interactions with GLU523 and LEU479. This study highlights CHEMBL3926946 as a robust lead for ROR1-targeted cancer therapy and emphasizes the utility of combining wt-MetaD and ABFEP for reliable hit prioritization. We employed a multilayered in silico pipeline integrating high-throughput virtual screening, long-timescale molecular dynamics (MD), wt-MetaD, and ABFEP. Ligands and protein were prepared using the OPLS2005 force field, and all stages up to wt-MetaD were conducted in Maestro (v12.8.117) using the same force field. A library of ~ 4 million compounds yielded 137 candidates, Further shortlisted via MD. 7 high-confidence molecules underwent 5 independent MD replicates with randomized seeds to ensure statistical robustness. The top 3 compounds were validated by 1 μs (1000 ns) simulations to assess long-term conformational stability and wt-MetaD to reveal deep minimum. ABFEP calculations were performed using the CGenFF force field in NAMD 3.0. We benchmarked ABFEP protocol against experimentally validated ligands, successfully reproducing experimental binding free energies (ΔG), confirming the protocol's predictive reliability.

Neural network potentials trained on quantum-mechanical data can calculate molecular interactions with relatively high speed and accuracy. However, not all neural network potentials are suitable for molecular simulations, as they might exhibit instabilities, nonphysical behavior, or lack accuracy. To assess the reliability of neural network potentials, a series of tests is conducted during model training, in the gas phase, and in the condensed phase. The testing procedure is performed for eight in-house neural network potentials based on the ANI-2x data set, using both the ANI-2x and MACE architectures. This consistent framework allows an evaluation of the effect of the model architecture on its performance. For comparison, we also perform stability tests of the publicly available neural network potentials: ANI-2x, ANI-1ccx, MACE-OFF23, and AIMNet2. The results show that the different models have different weaknesses. A normal-mode analysis of 14 simple benchmark molecules with large displacements from the energy minima revealed that the published MACE-OFF23-S model shows large deviations from the reference quantum-mechanical energy surface. Also, some MACE models with a reduced number of parameters failed to produce stable molecular dynamics simulations in the gas phase, and all MACE models exhibit unfavorable behavior during steric clashes. In addition, the published ANI-2x and one of the in-house MACE models are not able to reproduce the structure of liquid water at ambient conditions, forming an amorphous solid phase instead. For the ANI-1ccx model, the multibody interactions in the condensed water phase lead to nonphysical additional energy minima in bond length and bond angle space, which caused a phase transition to an amorphous solid. Out of all 13 considered public and in-house models, only one in-house model based on the ANI-2x B97-3c data set shows better agreement with the experimental radial distribution function of water than the simple molecular mechanics TIP3P and OPC models. Protein-ligand interaction energies for the four benchmark systems TYK2, CDK2, JNK1, and P38 show that almost all models exhibit a higher correlation with experimental binding affinities than the Chemgauss4 docking score (average R2 > 0.16). With an average R2 of 0.43, the ANI-2x model outperforms molecular mechanics calculations with the GAFF2 force field and DFTB3 semiempirical calculations (average R2 of 0.39 and 0.38), approaching the accuracy of absolute binding free energy calculations (average R2 of 0.52). However, the rather mixed results for the different machine learning potentials show that great care must be taken during model training and when selecting a neural network potential for real-world applications.

Benzimidazole hybrids are promising antibacterial agents, but the growing problem of antibiotic resistance has led to the necessity of developing novel compounds with enhanced antimicrobial activity. This study utilizes AI methods to generate new antibacterial compounds based on benzimidazole derivatives. We compiled a data set of these hybrids to explore their chemical space and identify effective scaffolds. An interpretable machine learning model was trained, achieving an R2 of 0.81 and RMSE of 0.212 for bioactivity prediction. Additionally, we employed a reinforcement learning model to create novel hybrid antibiotics through fragment-based, linker-based, and de novo approaches, selecting candidates based on bioactivity and off-target effects. This process yielded 56 novel synthetically feasible compounds with lower minimum inhibitory concentrations and improved drug-like properties. Molecular docking studies and absolute binding free energy (ABFE) calculations revealed that these generated molecules exhibit higher binding affinities to target proteins compared to approved antibiotics like ciprofloxacin and novobiocin.

Protein-ligand binding affinity prediction is a key aspect in computational small molecule drug discovery. Several recent studies have demonstrated that molecular simulations based on alchemical absolute binding free-energy (ABFE) calculations are an accurate and broadly applicable tool for this purpose. However, the use of current ABFE protocols in large scale drug discovery projects occasionally leads to unstable simulations and poor convergence. To address these problems, we have implemented several optimizations of the ABFE protocol. First, a new algorithm to choose the protein-ligand pose restraints was developed to prevent numerical instabilities. By considering protein-ligand hydrogen bonds, the restraint selection now incorporates data on the key interactions to improve the convergence. Second, an optimization of the annihilation protocol was conducted to minimize the resulting error of the free energy. Third, a rearrangement of the order with which interactions (electrostatics, Lennard-Jones, restraints, intramolecular torsions) are scaled resulted in a systematic improvement of the precision. The results from four protein-ligand benchmark systems (TYK2, P38, JNK1, and CDK2) show significantly lower variances of the free energy results and improvements of up to 0.23 kcal/mol for the root-mean-square error, compared to the original protocol.

Free-energy perturbation (FEP) calculations have emerged as a promising tool for the accurate prediction of ligand binding affinities. However, their widespread adoption in drug discovery pipelines has been hindered by long computation times and complex workflow setups. Here, we introduce an optimized graphics processing unit (GPU)-resident FEP implementation within GROMACS. The GPU-enabled FEP calculations are validated on a benchmark system containing eight ligand-protein pairs, including two charged ligands, on both the Nvidia A100 and the MetaX C500 GPU platforms. The absolute binding free energies predicted on these GPUs show excellent agreement (around 1.0 kcal/mol) with previous CPU-computed results. Compared to a 32-core CPU implementation, the GPU-accelerated FEP calculations demonstrate significant speed-ups, with up to nearly 800 and 400% improvements on Nvidia A100 and MetaX C500 GPUs, respectively. The end-to-end absolute binding free-energy calculations for the benchmark systems are reduced from 400 h to around 48 h on the A100 GPU. These advancements aim to provide the alchemical free-energy community with a fast and efficient way of conducting FEP calculations, thereby paving the way for a highly accurate and computationally efficient solution in predicting ligand-protein binding free energies. All codes, data, and scripts are included in our open-source project, FEP-on-GPU workflow, freely available at https://github.com/yiqichenshallwetalk/FEP-on-GPU-Workflow.

Identification of novel inhibitors for epidermal growth factor receptor tyrosine kinase using absolute binding free-energy simulations.

Apoptosis signal-regulating kinase 1 (ASK1), a key component of the mitogen-activated protein kinase (MAPK) cascades, has been identified as a promising therapeutic target owing to its critical role in signal transduction pathways. In this study, we proposed novel pyridin-2-yl urea inhibitors exhibiting favorable physicochemical properties. The potency of these compounds was validated through in vitro protein bioassays. The inhibition (IC50) of compound 2 was 1.55 ± 0.27 nM, which was comparable to the known clinical inhibitor, Selonsertib. To further optimize the hit compounds, two possible binding modes were initially predicted by molecular docking. Absolute binding free energy (BFE) calculations based on molecular dynamics simulations further discriminated the binding modes, presenting good tendency with bioassay results. This strategy, underpinned by BFE calculations, has the great potential to expedite the drug discovery process in the targeting of ASK1 kinase.

BPS2025 - Automated, efficient, and rigorous absolute binding free energy calculations

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.

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.

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