Water Networks in Complexes between Proteins and FDA-Approved Drugs

Water molecules at protein-ligand interfaces are often of significant pharmaceutical interest, owing in part to the entropy which can be released upon the displacement of an ordered water by a therapeutic compound. Protein structures may not, however, completely resolve all critical bound water molecules, or there may be no experimental data available. As such, predicting the location of water molecules in the absence of a crystal structure is important in the context of rational drug design. Grand canonical Monte Carlo (GCMC) is a computational technique that is gaining popularity for the simulation of buried water sites. In this work, we assess the ability of GCMC to accurately predict water binding locations, using a dataset that we have curated, containing 108 unique structures of complexes between proteins and Food and Drug Administration (FDA)-approved small-molecule drugs. We show that GCMC correctly predicts 81.4% of nonbulk crystallographic water sites to within 1.4 Å. However, our analysis demonstrates that the reported performance of water prediction methods is highly sensitive to the way in which the performance is measured. We also find that crystallographic water sites with more protein/ligand hydrogen bonds and stronger electron density are more reliably predicted by GCMC. An analysis of water networks revealed that more than half of the structures contain at least one ligand-contacting water network. In these cases, displacement of a water site by a ligand modification might yield unexpected results if the larger network is destabilized. Cooperative effects between waters should therefore be explicitly considered in structure-based drug design.

End-to-end differentiable construction of molecular mechanics force fields

Molecular mechanics (MM) potentials have long been a workhorse of computational chemistry. Leveraging accuracy and speed, these functional forms find use in a wide variety of applications in biomolecular modeling and drug discovery, from rapid virtual screening to detailed free energy calculations. Traditionally, MM potentials have relied on human-curated, inflexible, and poorly extensible discrete chemical perception rules (atom types) for applying parameters to small molecules or biopolymers, making it difficult to optimize both types and parameters to fit quantum chemical or physical property data. Here, we propose an alternative approach that uses graph neural networks to perceive chemical environments, producing continuous atom embeddings from which valence and nonbonded parameters can be predicted using invariance-preserving layers. Since all stages are built from smooth neural functions, the entire process-spanning chemical perception to parameter assignment-is modular and end-to-end differentiable with respect to model parameters, allowing new force fields to be easily constructed, extended, and applied to arbitrary molecules. We show that this approach is not only sufficiently expressive to reproduce legacy atom types, but that it can learn to accurately reproduce and extend existing molecular mechanics force fields. Trained with arbitrary loss functions, it can construct entirely new force fields self-consistently applicable to both biopolymers and small molecules directly from quantum chemical calculations, with superior fidelity than traditional atom or parameter typing schemes. When adapted to simultaneously fit partial charge models, espaloma delivers high-quality partial atomic charges orders of magnitude faster than current best-practices with low inaccuracy. When trained on the same quantum chemical small molecule dataset used to parameterize the Open Force Field ("Parsley") openff-1.2.0 small molecule force field augmented with a peptide dataset, the resulting espaloma model shows superior accuracy vis-á-vis experiments in computing relative alchemical free energy calculations for a popular benchmark. This approach is implemented in the free and open source package espaloma, available at https://github.com/choderalab/espaloma.

Best practices for constructing, preparing, and evaluating protein-ligand binding affinity benchmarks [Article v0.1]

Free energy calculations are rapidly becoming indispensable in structure-enabled drug discovery programs. As new methods, force fields, and implementations are developed, assessing their expected accuracy on real-world systems (benchmarking) becomes critical to provide users with an assessment of the accuracy expected when these methods are applied within their domain of applicability, and developers with a way to assess the expected impact of new methodologies. These assessments require construction of a benchmark-a set of well-prepared, high quality systems with corresponding experimental measurements designed to ensure the resulting calculations provide a realistic assessment of expected performance when these methods are deployed within their domains of applicability. To date, the community has not yet adopted a common standardized benchmark, and existing benchmark reports suffer from a myriad of issues, including poor data quality, limited statistical power, and statistically deficient analyses, all of which can conspire to produce benchmarks that are poorly predictive of real-world performance. Here, we address these issues by presenting guidelines for (1) curating experimental data to develop meaningful benchmark sets, (2) preparing benchmark inputs according to best practices to facilitate widespread adoption, and (3) analysis of the resulting predictions to enable statistically meaningful comparisons among methods and force fields. We highlight challenges and open questions that remain to be solved in these areas, as well as recommendations for the collection of new datasets that might optimally serve to measure progress as methods become systematically more reliable. Finally, we provide a curated, versioned, open, standardized benchmark set adherent to these standards (PLBenchmarks) and an open source toolkit for implementing standardized best practices assessments (arsenic) for the community to use as a standardized assessment tool. While our main focus is free energy methods based on molecular simulations, these guidelines should prove useful for assessment of the rapidly growing field of machine learning methods for affinity prediction as well.

Water molecules at protein-drug interfaces: computational prediction and analysis methods

The fundamental importance of water molecules at drug-protein interfaces is now widely recognised and a significant feature in structure-based drug design. Experimental methods for analysing the role of water in drug binding have many challenges, including the accurate location of bound water molecules in crystal structures, and problems in resolving specific water contributions to binding thermodynamics. Computational analyses of binding site water molecules provide an alternative, and in principle complete, structural and thermodynamic picture, and their use is now commonplace in the pharmaceutical industry. In this review, we describe the computational methodologies that are available and discuss their strengths and weaknesses. Additionally, we provide a critical analysis of the experimental data used to validate the methods, regarding the type and quality of experimental structural data. We also discuss some of the fundamental difficulties of each method and suggest directions for future study.

Rimantadine Binds to and Inhibits the Influenza A M2 Proton Channel without Enantiomeric Specificity

The influenza A M2 wild-type (WT) proton channel is the target of the anti-influenza drug rimantadine. Rimantadine has two enantiomers, though most investigations into drug binding and inhibition have used a racemic mixture. Solid-state NMR experiments using the full length-M2 WT have shown significant spectral differences that were interpreted to indicate tighter binding for (R)- vs (S)-rimantadine. However, it was unclear if this correlates with a functional difference in drug binding and inhibition. Using X-ray crystallography, we have determined that both (R)- and (S)-rimantadine bind to the M2 WT pore with slight differences in the hydration of each enantiomer. However, this does not result in a difference in potency or binding kinetics, as shown by similar values for kon, koff, and Kd in electrophysiological assays and for EC50 values in cellular assays. We concluded that the slight differences in hydration for the (R)- and (S)-rimantadine enantiomers are not relevant to drug binding or channel inhibition. To further explore the effect of the hydration of the M2 pore on binding affinity, the water structure was evaluated by grand canonical ensemble molecular dynamics simulations as a function of the chemical potential of the water. Initially, the two layers of ordered water molecules between the bound drug and the channel's gating His37 residues mask the drug's chirality. As the chemical potential becomes more unfavorable, the drug translocates down to the lower water layer, and the interaction becomes more sensitive to chirality. These studies suggest the feasibility of displacing the upper water layer and specifically recognizing the lower water layers in novel drugs.

Best Practices for Alchemical Free Energy Calculations [Article v1.0]

Alchemical free energy calculations are a useful tool for predicting free energy differences associated with the transfer of molecules from one environment to another. The hallmark of these methods is the use of "bridging" potential energy functions representing alchemical intermediate states that cannot exist as real chemical species. The data collected from these bridging alchemical thermodynamic states allows the efficient computation of transfer free energies (or differences in transfer free energies) with orders of magnitude less simulation time than simulating the transfer process directly. While these methods are highly flexible, care must be taken in avoiding common pitfalls to ensure that computed free energy differences can be robust and reproducible for the chosen force field, and that appropriate corrections are included to permit direct comparison with experimental data. In this paper, we review current best practices for several popular application domains of alchemical free energy calculations performed with equilibrium simulations, in particular relative and absolute small molecule binding free energy calculations to biomolecular targets.

Ligand Binding Free Energies with Adaptive Water Networks: Two-Dimensional Grand Canonical Alchemical Perturbations

Computational methods to calculate ligand binding affinities are increasing in popularity, due to improvements in simulation algorithms, computational resources, and easy-to-use software. However, issues can arise in relative ligand binding free energy simulations if the ligands considered have different active site water networks, as simulations are typically performed with a predetermined number of water molecules (fixed N ensembles) in preassigned locations. If an alchemical perturbation is attempted where the change should result in a different active site water network, the water molecules may not be able to adapt appropriately within the time scales of the simulations—particularly if the active site is occluded. By combining the grand canonical ensemble (μVT) to sample active site water molecules, with conventional alchemical free energy methods, the water network is able to dynamically adapt to the changing ligand. We refer to this approach as grand canonical alchemical perturbation (GCAP). In this work we demonstrate GCAP for two systems; Scytalone Dehydratase (SD) and Adenosine A_2A receptor. For both systems, GCAP is shown to perform well at reproducing experimental binding affinities. Calculating the relative binding affinities with a naïve, conventional attempt to solvate the active site illustrates how poor results can be if proper consideration of water molecules in occluded pockets is neglected. GCAP results are shown to be consistent with time-consuming double decoupling simulations. In addition, by obtaining the free energy surface for ligand perturbations, as a function of both the free energy coupling parameter and water chemical potential, it is possible to directly deconvolute the binding energetics in terms of protein–ligand direct interactions and protein binding site hydration.

Prediction of the Closed Conformation and Insights into the Mechanism of the Membrane Enzyme LpxR

Covalent modification of outer membrane lipids of Gram-negative bacteria can impact the ability of the bacterium to develop resistance to antibiotics as well as modulating the immune response of the host. The enzyme LpxR from Salmonella typhimurium is known to deacylate lipopolysaccharide molecules of the outer membrane; however, the mechanism of action is unknown. Here, we employ molecular dynamics and Monte Carlo simulations to study the conformational dynamics and substrate binding of LpxR in representative outer membrane models as well as detergent micelles. We examine the roles of conserved residues and provide an understanding of how LpxR binds its substrate. Our simulations predict that the catalytic H122 must be Nε-protonated for a single water molecule to occupy the space between it and the scissile bond, with a free binding energy of -8.5 kcal mol^-1. Furthermore, simulations of the protein within a micelle enable us to predict the structure of the putative "closed" protein. Our results highlight the need for including dynamics, a representative environment, and the consideration of multiple tautomeric and rotameric states of key residues in mechanistic studies; static structures alone do not tell the full story.

Replica-Exchange and Standard State Binding Free Energies with Grand Canonical Monte Carlo

The ability of grand canonical Monte Carlo (GCMC) to create and annihilate molecules in a given region greatly aids the identification of water sites and water binding free energies in protein cavities. However, acceptance rates without the application of biased moves can be low, resulting in large variations in the observed water occupancies. Here, we show that replica-exchange of the chemical potential significantly reduces the variance of the GCMC data. This improvement comes at a negligible increase in computational expense when simulations comprise of runs at different chemical potentials. Replica-exchange GCMC is also found to substantially increase the precision of water binding free energies as calculated with grand canonical integration, which has allowed us to address a missing standard state correction.

Polymorphism in crystals of bis(4-bromophenyl)fumaronitrile through vapour phase growth

We have linked polymorphic selectivity of crystals of bis(4-bromophenyl)fumaronitrile grown via physical vapour transport and their fluorescence to the lattice properties of each polymorph.