Benchmarking Force Field and the ANI Neural Network Potentials for the Torsional Potential Energy Surface of Biaryl Drug Fragments.

Many drug molecules contain biaryl fragments, resulting in a torsional barrier corresponding to rotation around the bond linking the aryls. The potential energy surfaces of these torsions vary significantly because of steric and electronic effects, ultimately affecting the relative stability of the molecular conformations in the protein-bound and solution states. Simulations of protein-ligand binding require accurate computational models to represent the intramolecular interactions to provide accurate predictions of the structure and dynamics of binding. In this article, we compare four force fields [generalized AMBER force field (GAFF), open force field (OpenFF), CHARMM general force field (CGenFF), optimized potentials for liquid simulations (OPLS)] and two neural network potentials (ANI-2x and ANI-1ccx) for their ability to predict the torsional potential energy surfaces of 88 biaryls extracted from drug fragments. The root mean square deviation (rmsd) over the full potential energy surface and the mean absolute deviation of the torsion rotational barrier height (MADB) relative to high-level ab initio reference data (CCSD(T1)*) were used as the measure of accuracy. Uncertainties in these metrics due to the composition of the data set were estimated using bootstrap analysis. In comparison to high-level ab initio data, ANI-1ccx was most accurate for predicting the barrier height (rmsd: 0.5 ± 0.0 kcal/mol, MADB: 0.8 ± 0.1 kcal/mol), followed closely by ANI-2x (rmsd: 0.5 ± 0.0 kcal/mol, MADB: 1.0 ± 0.2 kcal/mol), then CGenFF (rmsd: 0.8 ± 0.1 kcal/mol, MADB: 1.3 ± 0.1 kcal/mol) and OpenFF (rmsd: 0.7 ± 0.1 kcal/mol, MADB: 1.3 ± 0.1 kcal/mol), then GAFF (rmsd: 1.2 ± 0.2 kcal/mol, MADB: 2.6 ± 0.5 kcal/mol), and finally OPLS (rmsd: 3.6 ± 0.3 kcal/mol, MADB: 3.6 ± 0.3 kcal/mol). Significantly, the neural network potentials (NNPs) are systematically more accurate and more reliable than any of the force fields. As a practical example, the NNP/molecular mechanics method was used to simulate the isomerization of ozanimod, a drug used for multiple sclerosis. Multinanosecond molecular dynamics (MD) simulations in an explicit aqueous solvent were performed, as well as umbrella sampling and adaptive biasing force-enhanced sampling techniques. The rate constant for this isomerization calculated using transition state theory was 4.30 × 10-1 ns-1, which is consistent with direct MD simulations.

Flexible Fitting of Small Molecules into Electron Microscopy Maps Using Molecular Dynamics Simulations with Neural Network Potentials.

Despite significant advances in resolution, the potential for cryo-electron microscopy (EM) to be used in determining the structures of protein-drug complexes remains unrealized. Determination of accurate structures and coordination of bound ligands necessitates simultaneous fitting of the models into the density envelopes, exhaustive sampling of the ligand geometries, and, most importantly, concomitant rearrangements in the side chains to optimize the binding energy changes. In this article, we present a flexible-fitting pipeline where molecular dynamics flexible fitting (MDFF) is used to refine structures of protein-ligand complexes from 3 to 5 Å electron density data. Enhanced sampling is employed to explore the binding pocket rearrangements. To provide a model that can accurately describe the conformational dynamics of the chemically diverse set of small-molecule drugs inside MDFF, we use QM/MM and neural-network potential (NNP)/MM models of protein-ligand complexes, where the ligand is represented using the QM or NNP model, and the protein is represented using established molecular mechanical force fields (e.g., CHARMM). This pipeline offers structures commensurate to or better than recently submitted high-resolution cryo-EM or X-ray models, even when given medium to low-resolution data as input. The use of the NNPs makes the algorithm more robust to the choice of search models, offering a radius of convergence of 6.5 Å for ligand structure determination. The quality of the predicted structures was also judged by density functional theory calculations of ligand strain energy. This strain potential energy is found to systematically decrease with better fitting to density and improved ligand coordination, indicating correct binding interactions. A computationally inexpensive protocol for computing strain energy is reported as part of the model analysis protocol that monitors both the ligand fit as well as model quality.

Simulating protein-ligand binding with neural network potentials.

Drug molecules adopt a range of conformations both in solution and in their protein-bound state. The strain and reduced flexibility of bound drugs can partially counter the intermolecular interactions that drive protein-ligand binding. To make accurate computational predictions of drug binding affinities, computational chemists have attempted to develop efficient empirical models of these interactions, although these methods are not always reliable. Machine learning has allowed the development of highly-accurate neural-network potentials (NNPs), which are capable of predicting the stability of molecular conformations with accuracy comparable to state-of-the-art quantum chemical calculations but at a billionth of the computational cost. Here, we demonstrate that these methods can be used to represent the intramolecular forces of protein-bound drugs within molecular dynamics simulations. These simulations are shown to be capable of predicting the protein-ligand binding pose and conformational component of the absolute Gibbs energy of binding for a set of drug molecules. Notably, the conformational energy for anti-cancer drug erlotinib binding to its target was found to be considerably overestimated by a molecular mechanical model, while the NNP predicts a more moderate value. Although the ANI-1ccX NNP was not trained to describe ionic molecules, reasonable binding poses are predicted for charged ligands, but this method is not suitable for modeling charged ligands in solution.