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

Hannah E Bruce Macdonald,Gregory A Ross,Jonathan W Essex,Christopher Cave-Ayland

doi:10.1021/acs.jctc.8b00614

Abstract

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 A2A 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.

Highlights

Issues arise in relative protein−ligand binding free energy calculations in cases where water molecules become trapped in the protein binding site
Conventional alchemical perturbation methods do not always cope with this situation, in occluded pockets, where exchange with bulk water may be prevented within a feasible time scale due to kinetic barriers
Grand canonical Monte Carlo (GCMC) has been combined with Multistate BAR (MBAR) to achieve dynamic adaptation of water networks with relative protein−ligand binding free energy calculations

Summary

■ INTRODUCTION

It is understood that active site water molecules can have a large impact on the binding free energy of a protein−ligand complex.[1−3] In recent years, efforts have been made to rationalize these active site interactions, addressing such questions as where are water molecules located, what impact do they have on ligand binding, and how can knowledge of the water molecules help to design new, high affinity molecules[4] Computational methods of varying speed and accuracy exist to locate active site water molecules and predict their binding affinities.[5−15] it is still unclear as to whether a ligand should be designed to displace an active site water molecule for entropic gain and direct interaction between protein and ligand, or if a ligand’s interactions should be optimized to utilize water molecules through stabilizing bridging interactions. For 1D-GCAP simulations, as only λ is varied and B is constant at the equilibrium Beq value, the relative free energy of two ligands can be determined using classical free energy approaches: thermodynamic integration (TI), Bennett’s Acceptance Ratio (BAR),[41] or Multistate BAR (MBAR).[42] As with running GCMC at a single B value, 1D-GCAP is only able to determine the equilibrium number and location of water molecules, rather than the binding affinities of the water network. I is the index over all states, Ui is the potential energy according to the ith Hamiltonian, μi is the chemical potential of the ith state, and N is the occupancy of water molecules of state x This 2D-MBAR allows the free energy of the ligand perturbation to be calculated from the entire 2D-GCAP simulations, using statistically optimal contributions from all simulated states. The cavity near ring B is naıvely solvated using ProtoMS48 during the system set up

■ METHODS

■ CONCLUSION

Findings

■ REFERENCES