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.

Large-scale Virtual Screening (VS) campaigns of compound libraries can significantly speed up candidate selection in the early stages of drug discovery. The most promising drug candidates are identified by Scoring Functions (SFs), which enable VS campaigns to rank candidate compounds according to their estimated binding affinities. These SFs are typically trained on experimental data reflecting binding affinities (e.g., Dissociation Constant (Kd) values), commonly used as proxies for protein-ligand binding free energies. Because experimental reference data are often unavailable or collected using inconsistent techniques and/or procedures between laboratories, we developed two computational workflows that generate configurational ensembles of soluble protein-ligand complexes with Molecular Dynamics (MD) and compute the Absolute Binding Free Energies (ABFEs) of the sampled ligand binding poses with implicit-solvent calculations. The resulting consistent large-scale datasets of ABFEs address two complementary aspects of virtual screening: quantitative binding affinity estimation and binding pose assessment. Our Binding Affinity Prediction (BAP) workflow estimated protein-ligand binding affinities for 4000+ complexes from the PDBbind 2020 dataset. Our Pose Selector (PS) workflow computed non-convergence ABFEs from short Molecular Dynamics (MD) simulations, estimating the stability of 800,000+ related binding poses. To produce ABFE data at this scale, our free-energy workflows classify, check, and repair input structures of protein-ligand complexes in a fully automated fashion. The workflow scripts, molecular dynamics data, and ABFE labels are publicly available, creating an extendable database of reference values for the development of Scoring Functions for Large-Scale Virtual Screening campaigns.

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.

Accurate prediction of inhibitor selectivity across proteinparaloguesremains a central challenge in computational drug discovery. Here,we perform a comparative assessment of three computational methodsMolecularMechanics/Poisson–Boltzmann Surface Area (MM/PBSA), AbsoluteBinding Free Energy (ABFE) and Umbrella Sampling (US) calculationsintheir ability to recapitulate PARP1 versus PARP2 selectivity for eightclinically relevant PARP enzyme inhibitors used in ovarian, breast,and prostate tumors, among others. We demonstrate how MM/PBSA calculationsoffer rapid and qualitative insights but show pronounced sensitivityto the chosen static conformational pose, being particularly challengingfor ligands with subtle energetic differences between distinct proteinparalogues. In contrast, both ABFE and US calculations using atomisticmodels with explicit solvent result in substantially improved agreementwith experimental binding affinities. The ABFE method exhibits thestrongest quantitative correlation with experimental binding freeenergy differences, remarkably reproducing selectivity trends evenamong nearly isoenergetic complexes. Notably, our structural contactanalysis reveals how contact connectivity controls ligand selectivity,providing valuable mechanistic and molecular insight into the keyresidues that stabilize each inhibitor in both protein enzymes. Together,our multimethod computational study contributes to elucidating potentialchemical modifications across the ligand chemical space to enhancepotency and specificity, informing the future design and evaluationof selective inhibitors for precision oncology, including therapiestargeting homologous recombination-deficient cancers.

Computer simulations are increasingly significant in drug discovery, especially for predicting ligand affinities through free energy calculations. Absolute alchemical free energy calculations are vital for assessing ligand specificity and selectivity, aiding in differentiating between on-target efficacy and off-target effects. The main challenges in these calculations involve capturing conformational changes in proteins and ligands, predicting binding poses correctly, and modeling molecular interactions using well-parametrized force fields. Here, we evaluate how ab initio derived nonbonded force field parameters predict the specificity among nine BRD4 inhibitors and the selectivity of bromosporine across 22 bromodomains. We replaced nonbonded Open Force Field Sage 2.0.0 parameters with atomic charges, van der Waals radii, and dispersion coefficients obtained from Minimal Basis Iterative Stockholder (D-MBIS) atom partitioning of the polarized electron density, along with incorporating ligand polarization energies. Our ligand force field parameters demonstrated a mean unsigned error of 0.48 kcal/mol in predicting absolute binding free energy for nine BRD4 inhibitors, showing a strong correlation with experimental results. In the bromosporine selectivity set, predictive errors resulted mainly from docking-derived binding poses. By focusing solely on experimentally resolved apo- and holo structures, we accurately replicated experimental selectivity rankings for seven out of eight receptors, identifying those with the highest and lowest binding affinities.

Absolute Binding Free Energy (ABFE) methods are among the most accurate computational techniques for predicting protein-ligand binding affinities, but their utility is limited by the need for many simulations of alchemically modified intermediate states. We propose Direct Binding Free Energy (DBFE), an end-state ABFE method in implicit solvent that requires no alchemical intermediates. DBFE outperforms OBC2 double decoupling on a host-guest benchmark and performs comparably to OBC2 MM/GBSA on a protein-ligand benchmark. Since receptor and ligand simulations can be precomputed and amortized across compounds, DBFE requires only one complex simulation per ligand compared to the many lambda windows needed for double decoupling, making it a promising candidate for virtual screening workflows. We publicly release the code for this method at https://github.com/molecularmodelinglab/dbfe.

Hyperpolarization-activated cyclic nucleotide-gated (HCN) channels are are a family of voltage-gated, cyclic-nucleotide modulated Na+/K+ channels that regulate spontaneous rhythmic electrical activity in both the heart and the brain. Understanding differences in the responsiveness to cyclic adenosine monophosphate (cAMP) modulation between HCN isoforms would offer insight into the specific binding interactions that drive channel activation. Using all-atom molecular dynamics (MD) simulations and the free-energy perturbation (FEP) approach, we determined the absolute binding free energy of cAMP to the the cyclicnucleotide-binding domain (CNBD) of HCN isoforms 1–4. By studying the free-energy of ligand binding to the various isoforms of HCN, our study advances the understanding of HCN channel activation and modulation mechanisms. Overall, our work offers insight into explaining differences in channel sensitivity across the isoforms of HCN.

B-cell lymphoma 2 (Bcl-2) is a critical antiapoptotic protein and a prime therapeutic target in numerous cancers. The natural product gossypol is a known inhibitor of the Bcl-2 family, but the precise molecular details of its interaction remain elusive, hindering rational drug design efforts. In this study, we employed a comprehensive computational strategy, combining ensemble docking with advanced funnel metadynamics (FM) simulations, to elucidate the binding mechanism of gossypol to Bcl-2 at an atomic level. Our ensemble docking approach successfully predicted a consensus binding pose within the canonical BH3-mimetic groove. Subsequent FM simulations calculated an absolute binding free energy (ΔG) of -6.81 ± 0.86 kcal/mol, which shows reasonable quantitative agreement with the available experimental data. The reconstructed free-energy surface revealed a complex, multistep binding pathway involving a globally stable binding pose and several distinct, metastable intermediate states. Analysis of these states showed that hydrophobic forces are the primary drivers of binding. Furthermore, the interaction is markedly asymmetric; half of the gossypol molecule predominantly anchors the ligand into the P2 and P3 pockets in the most stable binding mode. Crucially, we demonstrate that gossypol binding reduces the overall flexibility of the binding site and that each binding state is characterized by a unique pattern of conformational stabilization across the four pockets. These findings provide an unprecedentedly detailed and dynamic roadmap of the gossypol-Bcl-2 interaction, offering crucial insights for the future structure-based design of next-generation inhibitors.

Free energy perturbation (FEP) is considered the gold-standard simulation method for estimating small molecule binding affinity, a quantity of vital importance to drug discovery. The accuracy of FEP critically depends on an accurate model of the protein-ligand complex as an initial condition for the underlying molecular dynamics simulation. This requirement has limited the impact of FEP in earlier stages of the discovery process, where appropriate experimental crystal structures are rarely available. The latest generation of structure prediction models, such as Boltz-2, promise to overcome this limitation by predicting protein-ligand complex structures. In this work, we combine Boltz-2 with our own absolute FEP protocol to build Boltz-ABFE, a robust pipeline for estimating the absolute binding free energies (ABFE) in the absence of experimental crystal structures. We investigate the quality of the structures predicted by Boltz-2, propose automated approaches to improve structures for use in molecular dynamics simulations, and demonstrate the effectiveness of the Boltz-ABFE pipeline for four protein targets from the FEP+ benchmark set. Demonstrating the feasibility of absolute FEP simulations without experimental crystal structures, Boltz-ABFE significantly expands the domain of applicability of FEP, paving the way toward accelerated early stage drug discovery via accurate, structure-based affinity estimation.

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.

ConspectusViral proteases are essential enzymes required for viral replication and assembly, making them prime antiviral drug targets. However, under the selective pressure of protease inhibitors, viruses can acquire mutations that reduce drug binding efficacy, posing significant challenges in both chronic infections (e.g., HIV, HCV) and acute infections like COVID-19, where mutations in the SARS-CoV-2 main protease (Mpro) have been reported to compromise the efficacy of drugs such as nirmatrelvir. To address these challenges, mainstream strategies in combating viral protease drug resistance mutations include combination therapies and targeting evolutionarily conserved regions of viral proteases. By disrupting multiple stages of the viral lifecycle or focusing on functionally indispensable residues, these strategies aim to develop next-generation antivirals that remain effective against evolving viral mutations.Here we provide an account of our laboratory's journey with our collaborators of the past 5 years, started during the COVID-19 pandemic and continued beyond it, in developing a multipronged strategy to combat antiviral drug resistance. The journey began with the discovery of compound 17 from screening our in-house α-ketoamide library with inhibitory activity against both the proteasome and SARS-CoV-2. We subsequently designed more specific covalent Mpro inhibitors with different warheads, such as H135 and H102 displaying potent activity, using computational and structural insights. H135 exhibited broad anti-SARS-CoV-2 activity, including alpha, delta, XBB.1.5, BA.5.2, EG.5.1, and JN.1.1 variants, in VeroE6-TMPRSS2 cells. Particularly, we observed an unusual distortion of the geometry of the catalytic dyad of Mpro in the X-ray crystal structure of H102, in which H102 induced conformational change of the catalytic residue His41. Using H102 as a model compound, we demonstrated that inducing conformational changes in His41 slightly enhances the antiresistance profile of inhibitors. Besides conventional protease inhibition, we attempted a new alternative strategy of protease degradation and developed HP211206, the first reported PROTAC molecule targeting SARS-CoV-2 Mpro, which is capable of degrading known drug-resistant Mpro mutants. To aid the design of these viral protease inhibitors and degraders, computational chemistry was used to develop an efficient method integrating the PDLD/S-LRA/β framework with quantum-mechanical calculations to evaluate both non-covalent and covalent contributions to absolute binding free energy between protein mutants and various inhibitors, which can help find inhibitors or degraders with high target binding. A vitality strategy was also employed to evaluate both the binding free energy of inhibitors (Ki) and the enzyme's catalytic efficiency parameters (kcat and KM) for natural substrates, which can help predict sites on the viral protease prone to drug resistance mutations and thus to be avoided in drug targeting. Additionally, a kinetic simulation framework was used to model the time dependence of inhibition (IC50(t)) for nirmatrelvir analogues, providing valuable insights into their time-dependent efficacy. Most recently, artificial intelligence was applied in the development of D2Screen, which incorporates deep learning into conventional virtual screening and enabled us to identify quinoline-based non-covalent Mpro inhibitors exhibiting anti-drug-resistance activity against the E166V Mpro mutant. Together, these synthetic, computational, structural, and biological studies illustrate a multipronged strategy for developing more effective therapeutics that are less susceptible to drug resistant mutation.

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.

Protein-protein (pp) interactions make up a varied class of potential therapeutic targets for malignancy, which cause disturbances in the alternating splicing of SR proteins by SRPK due to nuclear accumulation. SRPK1 phosphorylates several serine residues in the RS domain of ASF/SF2, a classic SR protein. Substrate selectivity is dependent on protein interactions beyond the kinase active site. The RS domain's phosphorylation cycle needs strong, sustained SRPK-SR binding. In light of this evidence, the present research investigates the dynamics of protein-protein complexes bound to nucleotide triphosphates and inhibitor molecules. The influence of protein-protein interactions and the role of small molecules at the binding site, exerting competitive inhibition, were examined using classical molecular dynamics simulations complemented by statistical analysis. Trajectory visualization and analysis revealed differential ASF-SRPK binding in the presence of ATP and MSC1186 (the most recent SRPK inhibitor). In fact, conformational changes and binding orientations of small molecules, together with ASF-SRPK interactions, are interdependent in sustaining biological functions. Free energy of binding of the small molecules at the active pocket was empirically obtained from the Generalized Born Implicit Solvation Model with details of the residue-wise contribution toward binding. A precise absolute binding free energy approach, based on streamlined alchemical free energy perturbation, was applied to evaluate the binding affinity at the competitive pocket and yielded results consistent with experimental data. Enhanced sampling method─"random accelerated molecular dynamics (RAMD)", including PMF evaluation and path analysis─was employed to explore the egression route and investigate the exit dynamics of small molecules from the binding pocket. This study lays the groundwork for novel therapeutic design methodologies while improving our molecular-level understanding of protein-protein-small-molecule interactions linked to carcinogenesis.

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.

Relative binding free energy (RBFE) calculations, widely used to predict the potencies of congeneric small molecules binding to a protein receptor, can greatly increase the efficiency of the hit-to-lead and lead optimization stages of the drug discovery process. Traditional RBFE methods, however, cannot be easily applied to small molecules lacking a common core or binding mode, precluding their use in a challenging but crucial component of many drug discovery campaigns. In principle, an absolute binding free energy (ABFE) method can be applied to such molecules, but ABFE often suffers from high computational cost and poor statistical convergence due to the large amount of additional sampling required when compared to RBFE. Here, we introduce core-hopping binding free energy (CBFE) calculations, a computationally efficient framework for the accurate determination of relative binding free energies between small molecules with different cores, leveraging several recently developed techniques such as Alchemical Enhanced Sampling (ACES) with optimized transformation pathways and flexible λ-spacing, as well as λ-dependent Boresch restraints. We benchmark the performance of CBFE across 4 protein systems consisting of 56 small molecules, and find that the results are consistent with RBFE for a congeneric series of ligands and offer considerable improvement in computational cost and precision relative to ABFE results for a series of small molecules with diverse cores and binding modes. All CBFE-related developments are fully implemented in the GPU-accelerated AMBER free energy module (pmemd.cuda) and are available as part of the latest official AMBER release.

Atom surface site Interaction Points (AIP) which werepreviouslyused to predict association constants for synthetic host–guestsystems has been extended to protein–ligand complexes. AIPdescriptions of protein binding sites were obtained by combining alibrary of precomputed AIP descriptors for all protein functionalgroups with a graph-based substructure matching algorithm. The correspondingAIP description of ligands was obtained directly by footprinting themolecular electrostatic potential surface calculated using densityfunctional theory. These AIP descriptions were projected onto X-raycrystal structures of protein–ligand complexes to identifypairs of AIPs that were sufficiently close in space to constitutean intermolecular interaction. The overall free energy of bindingwas calculated by summing the contributions of each AIP contact andassociated desolvation. Application to the 94 complexes involvinguncharged ligands in CASF benchmark data set showed that the methodachieves a Pearson correlation coefficient of 0.76 and an RMSD of11 kJ mol–1 for absolute free energies of binding.

Accurate prediction of protein-ligand binding free energies is critical yet computationally demanding in drug discovery. Alchemical free energy methods (AFEMs) offer high accuracy but suffer from significant computational costs and complex modeling setup, such as tuning the λ-schedule of alchemical transformation. While conventional deep learning (DL) models may instantly predict binding affinity, they often require a large training set and exhibit limited generalizability across chemical space. To address these challenges, we introduce LamNet, an alchemical-path-aware graph neural network. LamNet integrates endpoint molecular states and the bridging alchemical path (parametrized by λ) into a physics-informed representation learning framework, explicitly modeling free energy changes along a chosen thermodynamic transformation pathway. Trained on molecular-dynamics-simulated data along alchemical pathways and incorporating data reliability metrics, LamNet accurately predicts relative binding free energies and absolute binding free energies, and optimizes λ-schedules to improve traditional AFEM convergence. Evaluations on diverse datasets (463 ligands, 16 proteins) demonstrate that LamNet achieves superior or comparable performance to state-of-the-art methods, including traditional AFEM, but with up to 1000-fold acceleration. These findings establish LamNet as a generalizable, physics-grounded, and cost-effective tool that not only accelerates computations but also provides a novel framework for integrating rigorous computational physics into modern DL-driven drug discovery workflows.

Powered by dramatic advances in computer hardware, the advent of ultralarge make-on-demand virtual libraries, and a shift in small-molecule discovery toward more challenging targets with limited known actives, there has been a growing interest in the development of performant virtual screening methods that can reliably deliver novel hits. We report on a new method called Glide WS, that builds on our earlier efforts (WScore) to introduce an explicit representation of water structure and dynamics to an empirical scoring function suitable for high-throughput docking. This scoring function has been carefully tuned using absolute binding free energy perturbation calculations (ABFEP). Compared with Glide SP, Glide WS offers significant gains in the two primary tasks for molecular docking in drug discovery, pose prediction and virtual screening enrichment. For docking accuracy, Glide WS achieves a self-docking accuracy of 92% on a diverse set of 1477 protein ligand complexes as compared to 85% for Glide SP, using a criterion of 2.5 Å. We also demonstrate significantly improved virtual screening enrichment using a diverse data set covering of 38 targets together with three different computationally generated libraries of decoys, combined with standard known ChEMBL actives. We focus on ligands ranked in the top few percent of the database (the subset that is relevant to practical virtual screening efforts) and demonstrate that, along with improved enrichment of ChEMBL actives, Glide WS achieves a remarkable reduction in the number of poorly scoring decoys (as calibrated by ABFEP calculations), across a high percentage of targets, as compared to Glide SP. These results suggest that considerably higher hit rates will be observed, as compared to conventional rigid receptor docking, in practical virtual screening applications.

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.