Enhanced sampling methods, such as hybrid Monte Carlo/molecular dynamics (MC/MD), grand canonical MC/MD, nonequilibrium candidate MC/MD, etc., are widely used to sample the slow rearrangement of interfacial water molecules in binding free energy calculations. However, direct comparison of the accuracies of these methods is oftentimes challenging due to differences in studied systems, simulation parameters, and system setup across studies. To overcome this challenge, we introduce a novel and well-defined model system for benchmarking water sampling methods: a closed, custom-defined C 90 fullerene. Our C 90 fullerene inherently adopts a distorted oblate-spheroidal geometry rather than being spherical; essentially, C 90 is large enough to be almost a miniature, capped nanotube, unlike its C 60 "buckyball" sibling. Unlike conventional fullerenes with nonpolar, hydrophobic cavities, our custom-defined fullerene incorporates some level of polarity in the form of modest partial charges, creating a quadrupolar cavity that energetically favors water binding. In our quadrupolar fullerene, positive and negative partial charges are distributed around the polar and the equatorial regions, respectively, on its approximate oblate spheroidal geometry. Due to a significant free energy barrier imposed by the fullerene wall, water exchange between the cavity and bulk solvent is essentially impossible on MD time scales. To allow solvent water to equilibrate in the cavity, we introduce a solvent inlet, by performing Hamiltonian replica exchange (i.e., HREX) simulations, which reveal that the quadrupolar fullerene cavity can accommodate a maximum of two waters. From these HREX simulations, we also identify the water binding sites inside the cavity and their occupancies, which show that fullerene conformations with 2 or 1 waters inside the cavity are dominant, with nearly equal free energies. We further validate a nonequilibrium switching (NES) protocol for computing water displacement free energies, by removing a water from the conformations with 2 waters inside the cavity and comparing with the relevant free energies from the HREX simulations. Our NES free energies agree within statistical error with those obtained from HREX. Our findings establish the quadrupolar fullerene as a well-defined, reproducible test system for rigorously benchmarking current as well as future enhanced sampling techniques targeting slow water sampling in confined environments.

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.

Partial atomic charges are a fundamental component underlying classical molecular simulations, but assigning charges remains a computational bottleneck; many common methods rely on quantum mechanical calculations that scale poorly with molecular size and are sensitive to the choice of conformer geometry. We introduce Open Force Field (OpenFF) AshGC, a new graph convolutional neural network charge model, as well as the Sage 2.3.0 small molecule force field for drug-like molecules parametrized to be consistent with AshGC. AshGC is designed to efficiently produce conformer-independent charges of semiempirical quality at linear cost for molecules of all sizes, from small molecules to macromolecules. AshGC largely generates charges within the range of other accepted AM1-BCC backends such as OpenEye's oequacpac and AmberTools' sqm, deviating most in smaller molecules between 4 and 9 heavy atoms, in negatively charged molecules, and areas of chemistry underrepresented in the training set, such as particular sulfur- and phosphorus-containing functional groups. We further present the development and performance of Sage 2.3.0, which has both Lennard-Jones and valence parameters that are retrained to be consistent with neural network charges for the first time. Benchmarks spanning gas-phase geometry optimization through protein-ligand binding free energies show Sage 2.3.0 performs comparably to earlier Sage releases, with modest improvements in condensed-phase properties and a slight decrease in nonaqueous solvation free energy accuracy. As with other OpenFF force fields, Sage 2.3.0 was validated in protein-ligand benchmarks to be compatible with Amber protein force fields. All data are publicly available, along with scripts and environments for reproducing the training and benchmarking of AshGC and Sage 2.3.0 at https://github.com/openforcefield/ashgc-v1.0-fit and https://github.com/openforcefield/ash-sage-rc2 respectively.

Physics-based methods such as protein-ligand binding free energy calculations have been increasingly adopted in early-stage drug discovery to prioritize promising compounds for synthesis. However, the accuracy of these methods is highly dependent on details of the calculation and choices made while preparing the ligands and protein ahead of running calculations. During ligand preparation, researchers typically assign partial atomic charges to each ligand atom using a specific ligand conformation for charge assignment, often the input conformer. While it is a well known problem that partial charge assignment is dependent on conformation, little investigation has explored the downstream effects of varied partial charge assignment on free energy estimates. Preliminary benchmarks from Open Free Energy Project show that generating partial charges from different input conformers leads to variation of up to ±5.3 kcal/mol in calculated relative binding free energies due to variation in partial charges alone. In this study, we more systematically explore this issue, investigating systems where differences in partial charge generation (such as those caused by input conformer choice, partial charge engine, and hardware) may lead to differences in calculated absolute hydration free energy (AHFE) values. We demonstrate that supplying different input conformers to a partial charge engine can result in atomic partial charge discrepancies of up to 0.681 e, resulting in differences in calculated AHFE of 6.9 ± 0.1 kcal/mol. We find that even relatively small variations in partial charge assignment can result in notable differences in calculated AHFE, and thus care should be taken when assigning partial charges to ensure reproducibility and accuracy of any resulting free energy calculations.

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

Organophosphorus (OP) compounds are among the most toxic of chemical substances and widely used as insecticides, pesticides, and chemical warfare agents. The most important enzyme inhibited by OP compounds is acetylcholinesterase (AChe). Inactivation of AChe function results in the accumulation of neurotransmitter, leading to death due to serious respiratory disorders. Organophosphorus hydrolase (OPH), also called phosphotriesterase, is a homo-dimeric metalloenzyme that can hydrolyze various OP agents in the circulatory system, resulting in products that are generally of reduced toxicity. The best OPH substrate found to date is the insecticide diethyl p-nitrophenyl phosphate (paraoxon). Most structural and kinetic studies assume that the binding orientation of paraoxon is identical to that of diethyl 4-methylbenzylphosphonate, which is the only substrate analog co-crystallized with OPH. In the current work, we used a combined docking and molecular dynamics (MD) approach to predict the likely binding mode of paraoxon in the OPH active site. We identified a potential binding mode of paraoxon that does not match the binding mode of diethyl 4-methylbenzylphosphonate. Then, we used the predicted binding mode to run MD simulations on the wild type (WT) OPH complexed with paraoxon, and OPH mutants complexed with paraoxon. Additionally, we identified 3 hot-spot residues (D253, H254, and I255) involved in the stability of the OPH active site. To further assess these predictions, we then experimentally assayed single and double mutants involving these residues (D253E, H254S, I255S, D253E-H254R and D253E-I255G) for hydrolytic activity against paraoxon. Computational structural analysis of protein-substrate dynamics shows different hydrogen bonding profiles for mutants involving D253 (D253E, D253E-H254R, and D253E-I255G) compared to WT OPH. Additionally, the binding free energy calculations and the experimental kinetics (particularly, kcat and KM) of the reactions between each OPH mutant and paraoxon show that mutated forms D253E, D253E-H254R, and D253E-I255G exhibit enhanced activity over WT OPH. Interestingly, our experimental results show that the activity of the double mutant D253E-H254R increased by 19-fold compared to WT OPH.

Regulator of G protein signaling 2 (RGS2) negatively modulates signaling downstream of G protein-coupled receptors by accelerating GTP hydrolysis at Gα subunits of heterotrimeric G proteins. Decreased RGS2 levels are implicated in numerous diseases, including cardiovascular disease and asthma. Thus, identifying selective means of enhancing RGS2 protein levels would be a viable therapeutic strategy. RGS2 is rapidly degraded through the ubiquitin-proteasomal pathway, and we previously identified F-box only protein 44 (FBXO44) as the substrate recognition component of the E3 ligase responsible for facilitating RGS2 degradation. As such, the RGS2-FBXO44 interaction is a potential target for pharmacological intervention. Detailed information on the FBXO44 recognition site (degron) in RGS2 will aid in structure-based small-molecule inhibitor design, as well as in identifying additional FBXO44 targets, which would help predict possible side effects of targeting this interaction. Thus, the goal of this study was to dissect the molecular properties for FBXO44 binding of the RGS2 degron. We used a peptide array utilizing systematic residue substitution, combined with AlphaFold modeling and molecular dynamics simulations, to identify several amino acid changes that altered binding both positively and negatively. Finally, we experimentally confirmed our results in cells through coimmunoprecipitation and proteasomal inhibition, using full-length RGS2. Altogether, these results provide structural insights into RGS2-FBXO44 binding, which will aid in structure-guided drug discovery efforts. It also provides a framework for building a consensus recognition motif for FBXO44, which could aid in identifying more substrates for this understudied F-box protein.

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

This review article provides an overview of structurally oriented experimental datasets that can be used to benchmark protein force fields, focusing on data generated by nuclear magnetic resonance (NMR) spectroscopy and room temperature (RT) protein crystallography. We discuss what the observables are, what they tell us about structure and dynamics, what makes them useful for assessing force field accuracy, and how they can be connected to molecular dynamics simulations carried out using the force field one wishes to benchmark. We also touch on statistical issues that arise when comparing simulations with experiment. We hope this article will be particularly useful to computational researchers and trainees who develop, benchmark, or use protein force fields for molecular simulations.

The formation of protein-ligand complexes involves displacement of water molecules that were previously occupying the protein's binding site. In some cases, however, some water molecules may not be displaced by the ligand's binding, and they can stabilize the complex by mediating the interactions between the ligand and the protein. A relative binding free energy (RBFE) calculation between two ligands, one of which binds to the protein with an intermediate water while the other displaces the water, can yield wrong results if the water fails to rearrange itself within the simulation timescale. Enhanced sampling methods have previously been used to address the sampling of such "trapped" waters, inserting or deleting waters in the protein's binding site during ligand transformation. While sometimes effective, the enhanced sampling methods typically require long simulation times to converge and may lead to differences in RBFE estimates (i.e., hysteresis) based on initial water placement. In this study, we present a non-equilibrium switching (NES) method to calculate RBFEs in systems with trapped waters. Our approach requires the knowledge of the positions of the trapped waters prior to performing the free energy calculation for ligand transformation and then uses this information to efficiently calculate the RBFE between the ligands. In our simulation protocol, we perform ligand transformation in the binding site of the target protein by using three consecutive NES switches. The three NES switches implement restraints, transform the ligand, and then remove the restraints. We demonstrate that our NES simulation-based method results in RBFE estimates within 1.1 kcal mol-1 of experimental RBFEs, with associated statistical errors under 0.4 kcal mol-1, for eight systems involving trapped water displacement. Our method provides a computationally inexpensive alternative for estimating RBFEs for systems involving trapped waters by leveraging distributed computational resources.

We report the results of the SAMPL9 host-guest blind challenge for predicting binding free energies. The challenge focused on macrocycles from pillar[n]-arene and cyclodextrin host families, including WP6, and bCD and HbCD. A variety of methods were used by participants to submit binding free energy predictions. A machine learning approach based on molecular descriptors achieved the highest accuracy (RMSE of 2.04 kcal/mol) among ranked methods in the WP6 dataset. Interestingly, predictions for WP6 obtained via docking tended to outperform all methods (RMSE of 1.70 kcal/mol), most of which are MD based and computationally more expensive. In general, methods applying force fields achieved better correlation with experiments for WP6 opposed to the machine learning and docking models. In the cyclodextrin-phenothiazine challenge, the ATM approach emerged as the top performing method with RMSE less than 1.86 kcal/mol. Correlation metrics of ranked methods in this dataset was relatively poor compared to WP6. We also highlight several lessons learned to guide future work and help improve studies on the systems discussed. For example, WP6 may be present in other microstates other than its -12 state in the presence of certain guests. Machine learning approaches can be used to fine tune or help train force fields for certain chemistry (i.e WP6-G4). Certain phenothiazines occupy distinct primary and secondary orientations, some of which were considered individually for accurate binding free energies. The accuracy of predictions from certain methods while starting from a single binding pose/orientation demonstrate the sensitivity of calculated binding free energies to the orientation, and in some cases the likely dominant orientation for the system. Computational and experimental results suggest that guests phenothiazine core traverses both secondary and primary faces of the cyclodextrin hosts, bulky catioinic side chain will primarily occupy the primary face, and the phenothiazine core substituent resides at the larger secondary face.

Force fields are a key component of physics-based molecular modeling, describing the energies and forces in a molecular system as a function of the positions of the atoms and molecules involved. Here, we provide a review and scientific status report on the work of the Open Force Field (OpenFF) Initiative, which focuses on the science, infrastructure and data required to build the next generation of biomolecular force fields. We introduce the OpenFF Initiative and the related OpenFF Consortium, describe its approach to force field development and software, and discuss accomplishments to date as well as future plans. OpenFF releases both software and data under open and permissive licensing agreements to enable rapid application, validation, extension, and modification of its force fields and software tools. We discuss lessons learned to date in this new approach to force field development. We also highlight ways that other force field researchers can get involved, as well as some recent successes of outside researchers taking advantage of OpenFF tools and data.

We report the results of the SAMPL9 host–guest blind challenge for predicting binding free energies.

A wide range of density functional methods and basis sets are available to derive the electronic structure and properties of molecules. Quantum mechanical calculations are too computationally intensive for routine simulation of molecules in the condensed phase, prompting the development of computationally efficient force fields based on quantum mechanical data. Parametrizing general force fields, which cover a vast chemical space, necessitates generating sizable quantum mechanical datasets with optimized geometries and torsion scans. To achieve this efficiently, it is crucial to choose a quantum mechanical method that balances computational cost and accuracy. In this study we seek to assess the accuracy of quantum mechanical theory for specific properties such as conformer energies, electrostatic properties, and torsion energetics. To comprehensively evaluate various methods, we focus on a representative set of 59 diverse small molecules, comparing approximately 24 combinations of functional and basis sets against the reference level coupled cluster calculations at complete basis set limit.

We report the synthesis and characterization of sulfated pillar[5]arene hosts (P5S2-P5S10) that differ in the number of sulfate substituents. All five P5Sn hosts display high solubility in water (73-131 mM) and do not undergo significant self-association according to 1H NMR dilution experiments. The x-ray crystal structures of P5S6, P5S6 ⋅ Me6HDA, P5S8 ⋅ Me6HDA, and P5S10 ⋅ Me6HDA reveal one intracavity molecule of Me6HDA and several external molecules of Me6HDA which form a network of close methonium ⋅ ⋅ ⋅ sulfate interactions. The thermodynamic parameters of complexation between P5Sn and the panel of guests was measured by direct or competitive isothermal titration calorimetry. We find that the binding free energy toward a guest becomes more negative as the number of sulfate substituents increase. Conversely, the binding free energy of a specific sulfated pillar[5]arene toward a homologous series of guests becomes more negative as the number of NMe groups increases. The ability to tune the host ⋅ guest affinity by changing the number of sulfate substituents will be valuable in supramolecular polymers, separation materials, and latching applications.

As a model system, the binding pocket of the L99A mutant of T4 lysozyme has been the subject of numerous computational free energy studies. However, previous studies have failed to fully sample and account for the observed changes in the binding pocket of T4 L99A upon binding of a congeneric ligand series, limiting the accuracy of results. In this work, we resolve the closed, intermediate, and open states for T4 L99A previously reported in experiment in MD and establish definitions for these states based on the dynamics of the system. From this analysis, we arrive at two primary conclusions. First, assignment of simulation trajectories into discrete states should not be done simply based on RMSD to crystal structures as this can result in misassignment of states. Second, the different metastable conformations studied here need to be carefully treated, as we estimate the time scales for conformational interconversion to be on the order of 102 to 103 ns─far longer than time scales for typical binding calculations. We conclude with a discussion on the need to develop enhanced sampling methods to generally account for significant changes in protein conformation due to relatively small ligand perturbations.

Obtaining accurate binding free energies from in silico screens has been a longstanding goal for the computational chemistry community. However, accuracy and computational cost are at odds with one another, limiting the utility of methods that perform this type of calculation. Many methods achieve massive scale by explicitly or implicitly assuming that the target protein adopts a single structure, or undergoes limited fluctuations around that structure, to minimize computational cost. Others simulate each protein-ligand complex of interest, accepting lower throughput in exchange for better predictions of binding affinities. Here, we present the PopShift framework for accounting for the ensemble of structures a protein adopts and their relative probabilities. Protein degrees of freedom are enumerated once, and then arbitrarily many molecules can be screened against this ensemble. Specifically, we use Markov state models (MSMs) as a compressed representation of a protein's thermodynamic ensemble. We start with a ligand-free MSM and then calculate how addition of a ligand shifts the populations of each protein conformational state based on the strength of the interaction between that protein conformation and the ligand. In this work we use docking to estimate the affinity between a given protein structure and ligand, but any estimator of binding affinities could be used in the PopShift framework. We test PopShift on the classic benchmark pocket T4 Lysozyme L99A. We find that PopShift is more accurate than common strategies, such as docking to a single structure and traditional ensemble docking-producing results that compare favorably with alchemical binding free energy calculations in terms of RMSE but not correlation - and may have a more favorable computational cost profile in some applications. In addition to predicting binding free energies and ligand poses, PopShift also provides insight into how the probability of different protein structures is shifted upon addition of various concentrations of ligand, providing a platform for predicting affinities and allosteric effects of ligand binding. Therefore, we expect PopShift will be valuable for hit finding and for providing insight into phenomena like allostery.

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

DNA-encoded library technology grants access to nearly infinite opportunities to explore the chemical structure space for drug discovery. Successful navigation depends on the design and synthesis of libraries with appropriate physicochemical properties (PCPs) and structural diversity while aligning with practical considerations. To this end, we analyze combinatorial library design constraints including the number of chemistry cycles, bond construction strategies, and building block (BB) class selection in pursuit of ideal library designs. We compare two-cycle library designs (amino acid + carboxylic acid, primary amine + carboxylic acid) in the context of PCPs and chemical space coverage, given different BB selection strategies and constraints. We find that broad availability of amines and acids is essential for enabling the widest exploration of chemical space. Surprisingly, cost is not a driving factor, and virtually, the same chemical space can be explored with "budget" BBs.

Methods for calculating the relative binding free energy (RBFE) between ligands to a target protein are gaining importance in the structure-based drug discovery domain, especially as methodological advances and automation improve accuracy and ease of use. In an RBFE calculation, the difference between the binding affinities of two ligands to a protein is calculated by transforming one ligand into another, in the protein-ligand complex, and in solvent. Alchemical binding free energy calculations are often used for such ligand transformations. Such calculations are not without challenges, however; for example, it can be challenging to handle interfacial waters when these play a crucial role in mediating protein-ligand binding. In some cases, the exchange of the interfacial waters with solvent water might be very infrequent in the course of typical molecular simulations, and such interfacial waters can be considered trapped on the simulation timescale. In these cases, RBFE calculation between two ligands, where one ligand binds with a trapped water while the other ligand displaces it, can result in inaccuracies if the surrounding water structure is not sampled adequately for both ligands. So far, a popular choice for treating the trapped waters in RBFE calculations is to combine free energy calculations with enhanced sampling methods that insert/delete waters in the binding site. Despite recent developments in the enhanced sampling methods, they can result in hysteresis in the RBFE estimate, depending on whether the simulations were started with or without the trapped waters. In this study, we introduce an alternative method, separation of states, to calculate the RBFE between ligand pairs where the ligands bind to the protein with different numbers/positions of trapped waters. The separation of states approach treats the sampling of the trapped waters separately from the free energy calculation of the ligand transformation. In our method, a trapped water in protein's binding site is decoupled from the system first, and the cavity created by its decoupling is stabilized. We then grow a larger ligand into this cavity-- a ligand that is known to displace the trapped water. In this study, we show that our method results in precise and accurate estimates of RBFEs for ligand pairs involving the rearrangement of trapped water via RBFE calculations for five such ligand pairs. We have optimized our simulation protocol to be suited for large distributed computational resources and have automated our RBFE calculation workflow.