We present a set of Grand Challenges for predictive modeling in small molecule drug discovery, with the goal of defining, prioritizing, and quantifying the areas where computation can have transformative impact. Rather than offering another broad survey of methods, this paper articulates specific scientific and technical problems that limit progress today and proposes measurable criteria by which advances can be judged. Our objective is to align researchers, investors, and industry leaders around the challenges that matter most for advancing real drug discovery programs. This work seeks to broaden participation in drug discovery by providing a clear roadmap for contributors from the growing areas of artificial intelligence, machine learning, robotics, high-performance computing, and quantum computing, where tools are rapidly advancing but are often disconnected from the practical realities of medicinal chemistry and pharmacology. The insights presented here draw on extensive discussions with experienced drug hunters, computational method developers, and thought leaders across biotech, pharma, software, and venture capital, as well as lessons learned from our own drug discovery efforts. While there is substantial enthusiasm (particularly around AI) for revolutionizing drug discovery, this moment demands sharper problem definition. Without clearly articulated challenges and performance metrics, computational innovation risks optimizing for benchmarks rather than for translational impact. We therefore identify Grand Challenges across four domains: Chemistry, Structure, Energy, and Pharmacology. For each domain, we present a common framework: a well-defined challenge, the underlying physical principles, its relevance to drug discovery decision-making, the current state of the field, and quantitative metrics that define progress. By grounding computational ambition in concrete scientific problems, we aim to catalyze sustained, measurable advances rather than episodic waves of enthusiasm.

The increasing importance and predictive power of modern molecular modeling, driven by physics- and machine learning-based methods, necessitates a new collaborative architecture to replace the isolated, traditional model of software development. The traditional approach often led to redundant engineering effort, high costs, and opaque systems that limit reproducibility, independent scrutiny, and scientific independence. Additionally, it results in taxpayer-funded research being left siloed in commercial tools where it cannot have as much impact as if it were returned to the general public. This perspective advocates for permissively licensed open source software as a scientific and economic multiplier by reducing the duplication of effort, enabling scientific validation of modeling tools, and frictionless experimentation with new ideas. Coordinated, multi-project consortia, such as Open Force Field, Open Free Energy, OpenFold, and OpenADMET have formed to collaboratively build shared computational infrastructure and release all methods under permissive licenses. The success of these large-scale efforts requires organizational structures that extend beyond code. The Open Molecular Software Foundation (OMSF), a US nonprofit, serves as a domain-specific institutional home and fiscal sponsor. By providing governance, administrative infrastructure, and dedicated research software engineers, OMSF aligns incentives across academic and industrial stakeholders. This framework enables a synergistic ecosystem where projects interoperate to accelerate innovation, eliminate duplication, and ensure long-term software sustainability, thereby creating durable foundations that elevate the entire molecular modeling community.

Tumor Necrosis Factor (TNF) is a trimeric cytokine that exists in soluble (sTNF) and membrane-bound (mTNF) forms, both of which regulate immune responses through interactions with their cognate receptors. sTNF has been the subject of extensive biophysical investigations, leading to a mechanistic model in which a symmetric arrangement of the trimer promotes receptor signaling, whereas asymmetric conformations inhibit this function. In contrast, the structural and energetic landscape of mTNF remains largely under-explored. Here, we combined multi-microsecond unbiased molecular dynamics and Metadynamics-Adaptive Biasing Force simulations to characterize mTNF embedded in a physiologically relevant lipid bilayer. We show that in the absence of membrane engagement, the extracellular domain (ECD) dynamically samples multiple asymmetric conformations similar to those observed in sTNF. The association of the ECD with the membrane, mediated primarily by the basic residues R78–R82, R107, R108, R120, and R207, restricts conformational heterogeneity and stabilizes the symmetric state. By quantifying the energetic effects of ECD–membrane association, we demonstrate that the symmetric state in the membrane-bound ECD is stabilized over asymmetric conformations to a greater extent than in sTNF. This energetic effect, together with the spatial confinement imposed by the lipid bilayer, may explain the previously reported reduction in affinity of certain biologics for mTNF. In conclusion, our work elucidates how the membrane influences the structure, dynamics, and energetics of TNF. These mechanistic insights could guide future efforts to design mTNF-selective inhibitors that account for both membrane constraints and the energetic modulation of the ECD conformational landscape. ### Competing Interest Statement The authors have declared no competing interest.

AIMS:Choice of first-in-human dose has critical implications for the safety of Phase I participants as well as the likelihood of reaching the therapeutic dose range during escalation. In this analysis, we present a population concentration-response modelling approach for selecting the Phase I starting dose for a novel stimulator of interferon response cGAMP interactor 1 (STING) agonist, SNX281. METHODS:Given the immune agonist mechanism of SNX281, we opted to select the starting dose according to the minimum anticipated biological effect level (MABEL). To determine the MABEL concentration, we fitted a population concentration-response model to cytokine induction data from an ex vivo whole blood assay. We selected a whole blood assay to obviate the need for free fraction scaling for this highly protein-bound drug. We used the population concentration-response model to estimate the lower 10th percentile for the 10% maximal interferon-β response concentration of SNX281, which was chosen as the MABEL concentration. We translated the ex vivo MABEL concentration to a human MABEL dose using a human pharmacokinetic projection based on allometric scaling from preclinical species. RESULTS:The human dose-peak concentration relationship projection fell within 2-fold of the clinical result. After the application of a safety factor, the MABEL dose was applied in the clinic and did not demonstrate dose-limiting toxicities. CONCLUSIONS:Our novel population modelling-based MABEL strategy for first-in-human dose selection resulted in successful clinical translation of a small molecule STING agonist.

Generative modeling with artificial intelligence (GenAI) offers an emerging approach to discover novel, efficacious, and safe drugs by enabling the systematic exploration of chemical space and to design molecules that are synthesizable while also having desirable drug properties. However, despite rapid progress in other industries, GenAI has yet to demonstrate clear and consistent value in prospective drug discovery applications. In this Perspective, we argue that the ultimate goal of generative chemistry is not just to generate "new" or "interesting" molecules, but to generate "beautiful" molecules─those that are therapeutically aligned with the program objectives and bring value beyond traditional approaches. We focus on five essential considerations for the successful applications of GenAI for drug discovery (GADD): 1) chemical synthesizability (accounting for time/cost constraints); 2) favorable ADMET (absorption, distribution, metabolism, excretion, and toxicity) properties; 3) desirable target-specific binding to modulate the biological mechanism of interest; 4) the construction of appropriate multiparameter optimization (MPO) functions to drive the GenAI toward the project objectives; and 5) human feedback from experienced drug hunters. Interestingly, defining the beauty of a molecule in a drug discovery program is not always obvious, being context-dependent as data emerge and priorities shift, making the role of expert human input indispensable. While MPO frameworks using complex desirability functions or Pareto optimization can help operationalize multifaceted project objectives, they cannot yet fully capture the nuanced judgment of experienced drug hunters. Reinforcement learning with human feedback (RLHF) offers a path to guide the GenAI toward therapeutically aligned molecules, just as RLHF played a pivotal role in training large language models (LLMs) like ChatGPT, especially in aligning the model's behavior with human expectations. While not responsible for the model's base knowledge, RLHF is essential in shaping how the model responds. In addition to RLHF, future progress in GADD will depend on better property prediction models and explainable systems that provide insights to expert drug hunters. "Beauty is in the eyes of the beholder"─for drug discovery, beauty is judged by experienced drug hunters and clinical success.

In recent years, reinforcement learning (RL) has emerged as a valuable tool in drug design, offering the potential to propose and optimize molecules with desired properties. However, striking a balance between capabilities, flexibility, reliability, and efficiency remains challenging due to the complexity of advanced RL algorithms and the significant reliance on specialized code. In this work, we introduce ACEGEN, a comprehensive and streamlined toolkit tailored for generative drug design, built using TorchRL, a modern RL library that offers thoroughly tested reusable components. We validate ACEGEN by benchmarking against other published generative modeling algorithms and show comparable or improved performance. We also show examples of ACEGEN applied in multiple drug discovery case studies. ACEGEN is accessible at https://github.com/acellera/acegen-open and available for use under the MIT license.

The Structure and TOpology Replica Molecular Mechanics (STORMM) code is a next-generation molecular simulation engine and associated libraries optimized for performance on fast, vectorized central processor units and graphics processing units (GPUs) with independent memory and tens of thousands of threads. STORMM is built to run thousands of independent molecular mechanical calculations on a single GPU with novel implementations that tune numerical precision, mathematical operations, and scarce on-chip memory resources to optimize throughput. The libraries are built around accessible classes with detailed documentation, supporting fine-grained parallelism and algorithm development as well as copying or swapping groups of systems on and off of the GPU. A primary intention of the STORMM libraries is to provide developers of atomic simulation methods with access to a high-performance molecular mechanics engine with extensive facilities to prototype and develop bespoke tools aimed toward drug discovery applications. In its present state, STORMM delivers molecular dynamics simulations of small molecules and small proteins in implicit solvent with tens to hundreds of times the throughput of conventional codes. The engineering paradigm transforms two of the most memory bandwidth-intensive aspects of condensed-phase dynamics, particle–mesh mapping, and valence interactions, into compute-bound problems for several times the scalability of existing programs. Numerical methods for compressing and streamlining the information present in stored coordinates and lookup tables are also presented, delivering improved accuracy over methods implemented in other molecular dynamics engines. The open-source code is released under the MIT license.

Targeted protein degradation (TPD) is emerging as a promising therapeutic approach for cancer and other diseases, with an increasing number of programs demonstrating its efficacy in human clinical trials. One notable method for TPD is Proteolysis Targeting Chimeras (PROTACs, or heterobifunctional degraders) that selectively degrade a protein of interest (POI) through E3-ligase induced ubiquitination followed by proteasomal degradation. PROTACs utilize a warhead-linker-ligand architecture to bring the POI (bound to the warhead) and the E3 ligase (bound to the ligand) into close proximity. The resulting non-native protein-protein interactions (PPIs) formed between the POI and E3 ligase lead to the formation of a stable POI-degrader-ligase ternary complex, enhancing cooperativity for TPD. A significant challenge in PROTAC design is the time-consuming and resource-intensive screening of the degrader linkers to induce favorable non-native PPIs between POI and E3 ligase. In this work, we present a physics-based computational protocol to systematically predict non-canonical and metastable PPI interfaces between an E3 ligase and a given POI, aiding in the design of linkers to stabilize the PROTAC ternary complex and enhance degradation. In our protocol, we build the non-Markovian dynamic model using the Integrative Generalized Master Equation (IGME) method from approximately 1.5 millisecond all-atom molecular dynamics (MD) simulations of linker-less encounter complex, to systematically explore the inherent PPIs between the oncogene homologue (KRAS) protein and the von Hippel-Lindau (VHL) E3 ligase. Our IGME model successfully revealed six metastable states each containing a different PPI interface. We selected three of these metastable states containing promising PPIs for linker design. Our selection criterion included the thermodynamic and kinetic stabilities of these PPIs and the accessibility of the linker to the solvent-exposed sites on the warheads and the E3 ligand. One of our selected PPIs closely matches a recent co-crystal PPI interface structure induced by an experimentally designed PROTAC with potent degradation efficacy. We anticipate that our IGME approach has significant potential for widespread application in predicting metastable POI-ligase encounter complex interfaces that can enable subsequent rational design of novel PROTACs.

The therapeutic approach of targeted protein degradation (TPD) is gaining momentum due to its potentially superior effects compared with protein inhibition. Recent advancements in the biotech and pharmaceutical sectors have led to the development of compounds that are currently in human trials, with some showing promising clinical results. However, the use of computational tools in TPD is still limited, as it has distinct characteristics compared with traditional computational drug design methods. TPD involves creating a ternary structure (protein-degrader-ligase) responsible for the biological function, such as ubiquitination and subsequent proteasomal degradation, which depends on the spatial orientation of the protein of interest (POI) relative to E2-loaded ubiquitin. Modeling this structure necessitates a unique blend of tools initially developed for small molecules (e.g., docking) and biologics (e.g., protein-protein interaction modeling). Additionally, degrader molecules, particularly heterobifunctional degraders, are generally larger than conventional small molecule drugs, leading to challenges in determining drug-like properties like solubility and permeability. Furthermore, the catalytic nature of TPD makes occupancy-based modeling insufficient. TPD consists of multiple interconnected yet distinct steps, such as POI binding, E3 ligase binding, ternary structure interactions, ubiquitination, and degradation, along with traditional small molecule properties. A comprehensive set of tools is needed to address the dynamic nature of the induced proximity ternary complex and its implications for ubiquitination. In this Perspective, we discuss the current state of computational tools for TPD. We start by describing the series of steps involved in the degradation process and the experimental methods used to characterize them. Then, we delve into a detailed analysis of the computational tools employed in TPD. We also present an integrative approach that has proven successful for degrader design and its impact on project decisions. Finally, we examine the future prospects of computational methods in TPD and the areas with the greatest potential for impact.

The Alchemical Transfer Method (ATM) is herein validated against the relative binding free energies of a diverse set of protein-ligand complexes. We employed a streamlined setup workflow, a bespoke force field, and the AToM-OpenMM software to compute the relative binding free energies (RBFE) of the benchmark set prepared by Schindler and collaborators at Merck KGaA. This benchmark set includes examples of standard small R-group ligand modifications as well as more challenging scenarios, such as large R-group changes, scaffold hopping, formal charge changes, and charge-shifting transformations. The novel coordinate perturbation scheme and a dual-topology approach of ATM address some of the challenges of single-topology alchemical relative binding free energy methods. Specifically, ATM eliminates the need for splitting electrostatic and Lennard-Jones interactions, atom mapping, defining ligand regions, and post-corrections for charge-changing perturbations. Thus, ATM is simpler and more broadly applicable than conventional alchemical methods, especially for scaffold-hopping and charge-changing transformations. Here, we performed well over 500 relative binding free energy calculations for eight protein targets and found that ATM achieves accuracy comparable to existing state-of-the-art methods, albeit with larger statistical fluctuations. We discuss insights into specific strengths and weaknesses of the ATM method that will inform future deployments. This study confirms that ATM is applicable as a production tool for relative binding free energy (RBFE) predictions across a wide range of perturbation types within a unified, open-source framework.

Conformational samplingof complex biomolecules is an emergingfrontier in drug discovery. Advances in lab-based structural biologyand related computational approaches like AlphaFold have made greatstrides in obtaining static protein structures for biologically relevanttargets. However, biology is in constant motion, and many importantbiological processes rely on conformationally driven events. Conventionalmolecular dynamics (MD) simulations run on standard hardware are impracticalfor many drug design projects, where conformationally driven biologicalevents can take microseconds to milliseconds or longer. An alternativeapproach is to focus the search on a limited region of conformationalspace defined by a putative reaction coordinate (i.e., path collectivevariable). The search space is typically limited by applying restraints,which can be guided by insights about the underlying biological processof interest. The challenge is striking a balance between the degreeto which the system is constrained and still allowing for naturalmotions along the path. A plethora of restraints exist to limit thesize of conformational search space, although each has drawbacks whensimulating complex biological motions. In this work, we present athree-stage procedure to construct realistic path collective variables(PCVs) and introduce a new kind of barrier restraint that is particularlywell suited for complex conformationally driven biological events,such as allosteric modulations and conformational signaling. The PCVpresented here is all-atom (as opposed to C-alpha or backbone only)and is derived from all-atom MD trajectory frames. The new restraintrelies on a barrier function (specifically, the scaled reciprocalfunction), which we show is particularly beneficial in the contextof molecular dynamics, where near-hard-wall restraints are neededwith zero tolerance to restraint violation. We have implemented ourPCV and barrier restraint within a hybrid sampling framework thatcombines well-tempered metadynamics and extended-Lagrangian adaptivebiasing force (meta-eABF). We use three particular examples of highpharmaceutical interest to demonstrate the value of this approach:(1) sampling the distance from ubiquitin to a protein of interestwithin the supramolecular cullin-RING ligase complex, (2) stabilizingthe wild-type conformation of the oncogenic mutant JAK2-V617F pseudokinasedomain, and (3) inducing an activated state of the stimulator of interferongenes (STING) protein observed upon ligand binding. For examples 2and 3, we present statistical analysis of meta-eABF free energy estimatesand, for each case, code for reproducing this work.

In N-body applications, the efficient evaluation of range-limited forces depends on applying certain constraints, including a cut-off radius and force symmetry (Newton's Third Law). When computing the pair-wise forces in parallel, finding the optimal mapping of particles and computations to memories and processors is surprisingly challenging, but can result in greatly reduced data movement and computation. Despite FPGAs having a distinct compute model (BRAMs/network/pipelines) from CPUs and ASICs, mappings on FPGAs have not previously been studied in depth: it was thought that the half-shell method was preferred. In this work, we find that the Manhattan method is sur-prisingly compatible with FPGA hardware. With the cache overlapping technique proposed in this paper, the ultra-fine-grained data access demanded by the Manhattan method can be satisfied, despite the fact that the memory blocks on FPGAs appear to be insufficiently fine-grained. We further demonstrate that, compared to the traditional baseline half-shell method, approximately a half of the filters (preprocessors) can be removed without performance degradation. For communication, the amount of data transferred can be reduced by 40% - 75% in the most common multi-FPGA scenarios. Moreover, data transfers are almost perfectly balanced along all directions, and the optimization requires only minimal hardware resources. The practical consequence is that nearly 2 x to 4 x the workload can be handled without upgrading the network connections between FPGAs. This is a critical finding given the relatively limited bandwidth available in many common accelerator boards and the strong-scaling applications to which FPGA clusters are being applied.

Conformational sampling of complex biomolecules is an emerging frontier in drug discovery. Indeed, advances in lab-based structural biology and related computational approaches like AlphaFold have made great strides in obtaining static protein structures. However, biology is in constant motion and many important biological processes rely on conformationally-driven events. Unrestrained molecular dynamics (MD) simulations require that the simulated time be comparable to the real time of the biological processes of interest, rendering pure MD impractical for many drug design projects, where conformationally-driven biological events can take microseconds to milliseconds or longer. An alternative approach is to accelerate the sampling of specific motions by applying restraints, guided by insights about the underlying biological process of interest. A plethora of restraints exist to limit the size of conformational search space, although each has drawbacks when simulating complex biological motions. In this work, we introduce a new kind of restraint for molecular dynamics simulations (MD) that is particularly well suited for complex conformationallydriven biological events, such as protein-ligand binding, allosteric modulations, conformational signalling, and membrane permeability. The new restraint, which relies on a barrier function (the scaled reciprocal function) is particularly beneficial to MD, where hard-wall restraints are needed with zero tolerance to restraint violation. We have implemented this restraint within a hybrid sampling framework that combines metadynamics and extended-Lagrangian adaptive biasing force (meta-eABF). We use two particular examples to demonstrate the value of this approach: (1) quantification of the approach of E3-loaded ubiquitin to a protein of interest as part of the Cullin ring ligase and (2) membrane permeability of heterobi-functional degrader molecules with a large degree of conformational flexibility. Future work will involve extension to additional systems and benchmarking of this approach compared with other methods.

Source code and structural PDBs that went into the Prediction of ternary complex structures using molecular dynamics and proteomics paper.

The protein STING (stimulator of interferon genes) is a central regulator of the innate immune system and plays an important role in antitumor immunity by inducing the production of cytokines such as type I interferon (IFN). Activation of STING stems from the selective recognition of endogenous cyclic dinucleotides (CDNs) by the large, polar, and flexible binding site, thus posing challenges to the design of small molecule agonists with drug-like physicochemical properties. In this work we present the design of SNX281, a small molecule STING agonist that functions through a unique self-dimerizing mechanism in the STING binding site, where the ligand dimer approximates the size and shape of a cyclic dinucleotide while maintaining drug-like small molecule properties. SNX281 exhibits systemic exposure, STING-mediated cytokine release, strong induction of type I IFN, potent in vivo antitumor activity, durable immune memory, and single-dose tumor elimination in mouse models via a C max -driven pharmacologic response. Bespoke computational methods – a combination of quantum mechanics, molecular dynamics, binding free energy simulations, and artificial intelligence – were developed during the course of the project to design SNX281 by explicitly accounting for the unique self-dimerization mechanism and the large-scale conformational change of the STING protein upon activation. Over the course of the project, we explored millions of virtual molecules while synthesizing and testing only 208 molecules in the lab. This work highlights the value of a multifaceted computationally-driven approach anchored by methods tailored to address target-specific problems encountered along the project progression from initial hit to the clinic.

Targeted protein degradation (TPD) is a promising approach in drug discovery for degrading proteins implicated in diseases. A key step in this process is the formation of a ternary complex where a heterobifunctional molecule induces proximity of an E3 ligase to a protein of interest (POI), thus facilitating ubiquitin transfer to the POI. In this work, we characterize 3 steps in the TPD process. (1) We simulate the ternary complex formation of SMARCA2 bromodomain and VHL E3 ligase by combining hydrogen-deuterium exchange mass spectrometry with weighted ensemble molecular dynamics (MD). (2) We characterize the conformational heterogeneity of the ternary complex using Hamiltonian replica exchange simulations and small-angle X-ray scattering. (3) We assess the ubiquitination of the POI in the context of the full Cullin-RING Ligase, confirming experimental ubiquitinomics results. Differences in degradation efficiency can be explained by the proximity of lysine residues on the POI relative to ubiquitin. The formation of ternary degrader-protein complexes is a key step in the targeted degradation of proteins of interest. Here, the authors explore the structure and dynamics of such complexes applying high-performance computer simulations augmented with experimental data.

To help cells cope with protein misfolding and aggregation, Hsp70 molecular chaperones selectively bind a variety of sequences ("selective promiscuity"). Statistical analyses from substrate-derived peptide arrays reveal that DnaK, the E. coli Hsp70, binds to sequences containing three to five branched hydrophobic residues, although otherwise the specific amino acids can vary considerably. Several high-resolution structures of the substrate -binding domain (SBD) of DnaK bound to peptides reveal a highly conserved configuration of the bound substrate and further suggest that the substrate-binding cleft consists of five largely independent sites for interaction with five consecutive substrate residues. Importantly, both substrate backbone orientations (N- to C- and C- to N-) allow essentially the same backbone hydrogen-bonding and side-chain interactions with the chaperone. In order to rationalize these observations, we performed atomistic molecular dynamics simulations to sample the interactions of all 20 amino acid side chains in each of the five sites of the chaperone in the context of the conserved substrate backbone configurations. The resulting interaction energetics provide the basis set for deriving a predictive model that we call Paladin (Physics-based model of DnaK-Substrate Binding). Trained using available peptide array data, Paladin can distinguish binders and nonbinders of DnaK with accuracy comparable to existing predictors and further predicts the detailed configuration of the bound sequence. Tested using existing DnaK-peptide structures, Paladin correctly predicted the binding register in 10 out of 13 substrate sequences that bind in the N- to C- orientation, and the binding orientation in 16 out of 22 sequences. The physical basis of the Paladin model provides insight into the origins of how Hsp70s bind substrates with a balance of selectivity and promiscuity. The approach described here can be extended to other Hsp70s where extensive peptide array data is not available.

With the current pandemic, the central role that Molecular Dynamics simulation (MD) plays in drug discovery makes advances in MD performance urgent. Recent work has demonstrated that among COTS devices only FPGA-centric clusters can scale beyond a few processors for relevant targets; other work has shown that single FPGA performance compares favorably to that of a GPU. In this study we demonstrate that an additional factor of 4× performance can be achieved which results in a factor of 5× speed up over a GPU. The problem addressed is that the designs of the last decade no longer scale when the number of processing pipelines grows from around ten to the hundreds. We begin by systematically evaluating existing work, exposing its flaws, and proposing a series of new design solutions. There are four major contributions. First, we address the massive routing problem by augmenting the design with three minimal networks in logic and latency. Second, we have developed a novel asynchronous out-of-order communication mechanism that removes nearly all bubbles from the routing networks. Third, we find that inverting the standard particle access algorithm results in improved locality and performance. Finally, we have created a custom numerical format that increases precision while saving space and logic.

Molecular Dynamics (MD) simulations play a central role in physics-driven drug discovery. MD applications often use the Particle Mesh Ewald (PME) algorithm to accelerate electrostatic force computations, but efficient parallelization has proven difficult due to the high communication requirements of distributed 3D FFTs. In this paper, we present the design and implementation of a scalable PME algorithm that runs on a cluster of Intel Stratix 10 FPGAs and can handle FFT sizes appropriate to address real-world drug discovery projects (grids up to 1283). To our knowledge, this is the first work to fully integrate all aspects of the PME algorithm (charge spreading, 3D FFT/IFFT, and force interpolation) within a distributed FPGA framework. The design is fully implemented with OpenCL for flexibility and ease of development and uses 100 Gbps links for direct FPGA-to-FPGA communications without the need for host interaction. We present experimental data up to 4 FPGAs (e.g., 206 microseconds per timestep for a 65536 atom simulation and 643 3D FFT), outperforming GPUs. Additionally, we discuss design scalability on clusters with differing topologies up to 64 FPGAs (with expected performance greater than all known GPU implementations) and integration with other hardware components to form a complete molecular dynamics application. We predict best-case performance of 6.6 microseconds per timestep on 64 FPGAs.

Alchemical binding free energy (BFE) calculations offer an efficient and thermodynamically rigorous approach to in silico binding affinity predictions. As a result of decades of methodological improvements and recent advances in computer technology, alchemical BFE calculations are now widely used in drug discovery research. They help guide the prioritization of candidate drug molecules by predicting their binding affinities for a biomolecular target of interest (and potentially selectivity against undesirable antitargets). Statistical variance associated with such calculations, however, may undermine the reliability of their predictions, introducing uncertainty both in ranking candidate molecules and in benchmarking their predictive accuracy. Here, we present a computational method that substantially improves the statistical precision in BFE calculations for a set of ligands binding to a common receptor by dynamically allocating computational resources to different BFE calculations according to an optimality objective established in a previous work from our group and extended in this work. Our method, termed Network Binding Free Energy (NetBFE), performs adaptive BFE calculations in iterations, re-optimizing the allocations in each iteration based on the statistical variances estimated from previous iterations. Using examples of NetBFE calculations for protein binding of congeneric ligand series, we demonstrate that NetBFE approaches the optimal allocation in a small number (≤5) of iterations and that NetBFE reduces the statistical variance in the BFE estimates by approximately a factor of 2 when compared to a previously published and widely used allocation method at the same total computational cost.