Collaborative Assessment of Molecular Geometries and Energies from the Open Force Field

Machine-learned molecular mechanics force fields from large-scale quantum chemical data

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.

The Open Force Field Initiative: Open Software and Open Science for Molecular Modeling

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.

Complexation Mechanisms of Aqueous Amylose: Molecular Dynamics Study Using 3-Pentadecylphenol

Molecular dynamics (MD) simulations of linear amylose fragments containing 10 to 40 glucose units were used to study the complexation of the prototypical compound, 3-pentadecylphenol (PDP)─a natural product with surfactant-like properties─in aqueous solution. The amylose-PDP binding leverages mainly hydrophobic interactions together with excluded volume effects. It was found that while the most stable complexes contained PDP inside the helical structure of the amylose in the expected guest-host (inclusion) complexation manner, at higher temperatures, the commonly observed PDP-amylose complexes often involved more nonspecific interactions than inclusion complexation. In the case where a stoichiometric excess of PDP was added to the simulation box, self-aggregation of the small molecule precluded its ability to enter the internal helical part of the oligosaccharide, and as a result, inclusion complexation became ineffective. MD simulation trajectories were analyzed preliminarily using cluster analysis (CA), followed by more rigorous solvent accessible surface area (SASA) determination over the temperature range spanning from 277 to 433 K. It was found that using the SASA of PDP corrected for its intrinsic conformational changes, together with a generic hidden Markov model (HMM), an adequate quantification of the different types of PDP-amylose aggregates was obtained to allow further analysis. The enthalpy change associated with the guest-host binding equilibrium constant (Kgh) in aqueous solution was estimated to be -75 kJ/mol, which is about twice as high as one might expect based on experimentally measured values of similar complexes in the solid state where the (unsolvated) helical structure of amylose remains rigid. On the other hand, the nonspecific binding (Kns) enthalpy change associated with PDP-amylose interactions in the same solution environment was found to be about half of the inclusion complexation value.

Development and Comprehensive Benchmark of a High-Quality AMBER-Consistent Small Molecule Force Field with Broad Chemical Space Coverage for Molecular Modeling and Free Energy Calculation

Biomolecular simulations have become an essential tool in contemporary drug discovery, and molecular mechanics force fields (FFs) constitute its cornerstone. Developing a high quality and broad coverage general FF is a significant undertaking that requires substantial expert knowledge and computing resources, which is beyond the scope of general practitioners. Existing FFs originate from only a limited number of groups and organizations, and they either suffer from limited numbers of training sets, lower than desired quality because of oversimplified representations, or are costly for the molecular modeling community to access. To address these issues, in this work, we developed an AMBER-consistent small molecule FF with extensive chemical space coverage, and we provide Open Access parameters for the entire modeling community. To validate our FF, we carried out benchmarks of quantum mechanics (QM)/molecular mechanics conformer comparison and free energy perturbation calculations on several benchmark data sets. Our FF achieves a higher level of performance at reproducing QM energies and geometries than two popular open-source FFs, OpenFF2 and GAFF2. In relative binding free energy calculations for 31 protein-ligand data sets, comprising 1079 pairs of ligands, the new FF achieves an overall root-mean-square error of 1.19 kcal/mol for ΔΔG and 0.92 kcal/mol for ΔG on a subset of 463 ligands without bespoke fitting to the data sets. The results are on par with those of the leading commercial series of OPLS FFs.

Distribution of Bound Conformations in Conformational Ensembles for X-ray Ligands Predicted by the ANI-2X Machine Learning Potential

In this study, we systematically studied the energy distribution of bioactive conformations of small molecular ligands in their conformational ensembles using ANI-2X, a machine learning potential, in conjunction with one of our recently developed geometry optimization algorithms, known as a conjugate gradient with backtracking line search (CG-BS). We first evaluated the combination of these methods (ANI-2X/CG-BS) using two molecule sets. For the 231-molecule set, ab initio calculations were performed at both the ωB97X/6-31G(d) and B3LYP-D3BJ/DZVP levels for accuracy comparison, while for the 8,992-molecule set, ab initio calculations were carried out at the B3LYP-D3BJ/DZVP level. For each molecule in the two molecular sets, up to 10 conformations were generated, which diminish the influence of individual outliers on the performance evaluation. Encouraged by the performance of ANI-2x/CG-BS in these evaluations, we calculated the energy distributions using ANI-2x/CG-BS for more than 27,000 ligands in the protein data bank (PDB). Each ligand has at least one conformation bound to a biological molecule, and this ligand conformation is labeled as a bound conformation. Besides the bound conformations, up to 200 conformations were generated using OpenEye's Omega2 software (https://docs.eyesopen.com/applications/ omega/) for each conformation. We performed a statistical analysis of how the bound conformation energies are distributed in the ensembles for 17,197 PDB ligands that have their bound conformation energies within the energy ranges of the Omega2-generated conformation ensembles. We found that half of the ligands have their relative conformation energy lower than 2.91 kcal/mol for the bound conformations in comparison with the global conformations, and about 90% of the bound conformations are within 10 kcal/mol above the global conformation energies. This information is useful to guide the construction of libraries for shape-based virtual screening and to improve the docking algorithm to efficiently sample bound conformations.

Drug Aggregation of Sparingly-Soluble Ionizable Drugs: Molecular Dynamics Simulations of Papaverine and Prostaglandin F2α

Aggregation in aqueous solution can have important implications on both the in vivo exposure of a drug and its pharmaceutical manufacturability. However, the drug aggregates formed can be very small and, thus, difficult to interrogate experimentally. On the other hand, at higher supersaturations where larger aggregates are supported, the chemical system is inherently metastable and therefore likewise challenging to study from an experimental standpoint. Understanding aggregation behavior is further complicated in the case of ionizable drugs where, unlike neutral compounds, there can be uncertainty in the kinds of drug molecules (i.e., charged, neutral, or both) that become incorporated into various clusters, particularly at pH values near the pKa. In this paper, we apply physics-based all-atom molecular dynamics (MD) simulations to study aggregation in the weakly basic drug papaverine and in the weakly acidic drug prostaglandin F2α. We employ in silico tools to construct simulation workflows and comprehensive cluster analysis protocols to elucidate the size distributions and dynamics of the drug aggregates formed at both an experimentally relevant concentration and at high supersaturation. We build on a previously published treatment [Solubility of sparingly soluble ionizable drugs. Adv. Drug Deliv. Rev. 2007, 59, 568-590, DOI: 10.1016/j.addr.2007.05.008] to translate the predicted aggregate distributions of each ionized drug into corresponding pH-solubility curves that can be compared directly to experiment. Our findings show that the assumption of a single predominant (charged) aggregate can be misleading in interpreting experimental pH-solubility curves, as it does not adequately reflect the rich diversity revealed in our simulations. Beyond not accounting for the distribution of ionized drug-containing clusters actually observed in solution, for both drugs we find evidence that neutral drug molecules can also participate in the aggregation phenomena. Notably, we observe that many drug molecules remain as free monomers in solution even under simulated conditions designed to mimic those where there is significant deviation of the experimental pH-solubility curve from the Henderson-Hasselbalch (H-H) equation, often taken to be a clear signpost of drug aggregation.

Perspectives on Ligand/Protein Binding Kinetics Simulations: Force Fields, Machine Learning, Sampling, and User-Friendliness

Computational techniques applied to drug discovery have gained considerable popularity for their ability to filter potentially active drugs from inactive ones, reducing the time scale and costs of preclinical investigations. The main focus of these studies has historically been the search for compounds endowed with high affinity for a specific molecular target to ensure the formation of stable and long-lasting complexes. Recent evidence has also correlated the in vivo drug efficacy with its binding kinetics, thus opening new fascinating scenarios for ligand/protein binding kinetic simulations in drug discovery. The present article examines the state of the art in the field, providing a brief summary of the most popular and advanced ligand/protein binding kinetics techniques and evaluating their current limitations and the potential solutions to reach more accurate kinetic models. Particular emphasis is put on the need for a paradigm change in the present methodologies toward ligand and protein parametrization, the force field problem, characterization of the transition states, the sampling issue, and algorithms' performance, user-friendliness, and data openness.

A transferable double exponential potential for condensed phase simulations of small molecules

The Lennard-Jones potential is the most widely-used function for the description of non-bonded interactions in transferable force fields for the condensed phase. This is not because it has an optimal functional form, but rather it is a legacy resulting from when computational expense was a major consideration and this potential was particularly convenient numerically. At present, it persists because the effort that would be required to re-write molecular modelling software and train new force fields has, until now, been prohibitive. Here, we present Smirnoff-plugins as a flexible framework to extend the Open Force Field software stack to allow custom force field functional forms. We deploy Smirnoff-plugins with the automated Open Force Field infrastructure to train a transferable, small molecule force field based on the recently-proposed double exponential functional form, on over 1000 experimental condensed phase properties. Extensive testing of the resulting force field shows improvements in transfer free energies, with acceptable conformational energetics, run times and convergence properties compared to state-of-the-art Lennard-Jones based force fields.

Development and Benchmarking of Open Force Field 2.0.0: The Sage Small Molecule Force Field

We introduce the Open Force Field (OpenFF) 2.0.0 small molecule force field for drug-like molecules, code-named Sage, which builds upon our previous iteration, Parsley. OpenFF force fields are based on direct chemical perception, which generalizes easily to highly diverse sets of chemistries based on substructure queries. Like the previous OpenFF iterations, the Sage generation of OpenFF force fields was validated in protein-ligand simulations to be compatible with AMBER biopolymer force fields. In this work, we detail the methodology used to develop this force field, as well as the innovations and improvements introduced since the release of Parsley 1.0.0. One particularly significant feature of Sage is a set of improved Lennard-Jones (LJ) parameters retrained against condensed phase mixture data, the first refit of LJ parameters in the OpenFF small molecule force field line. Sage also includes valence parameters refit to a larger database of quantum chemical calculations than previous versions, as well as improvements in how this fitting is performed. Force field benchmarks show improvements in general metrics of performance against quantum chemistry reference data such as root-mean-square deviations (RMSD) of optimized conformer geometries, torsion fingerprint deviations (TFD), and improved relative conformer energetics (ΔΔE). We present a variety of benchmarks for these metrics against our previous force fields as well as in some cases other small molecule force fields. Sage also demonstrates improved performance in estimating physical properties, including comparison against experimental data from various thermodynamic databases for small molecule properties such as ΔHmix, ρ(x), ΔGsolv, and ΔGtrans. Additionally, we benchmarked against protein-ligand binding free energies (ΔGbind), where Sage yields results statistically similar to previous force fields. All the data is made publicly available along with complete details on how to reproduce the training results at https://github.com/openforcefield/openff-sage.