Predicting ligand binding affinity using on- and off-rates for the SAMPL6 SAMPLing challenge

Abstract

Interest in ligand binding kinetics has been growing rapidly, as it is being discovered in more and more systems that ligand residence time is the crucial factor governing drug efficacy. Many enhanced sampling methods have been developed with the goal of predicting ligand binding rates ([Formula: see text]) and/or ligand unbinding rates ([Formula: see text]) through explicit simulation of ligand binding pathways, and these methods work by very different mechanisms. Although there is not yet a blind challenge for ligand binding kinetics, here we take advantage of experimental measurements and rigorously computed benchmarks to compare estimates of [Formula: see text] calculated as the ratio of two rates: [Formula: see text]. These rates were determined using a new enhanced sampling method based on the weighted ensemble framework that we call "REVO": Reweighting of Ensembles by Variance Optimization. This is a further development of the WExplore enhanced sampling method, in which trajectory cloning and merging steps are guided not by the definition of sampling regions, but by maximizing trajectory variance. Here we obtain estimates of [Formula: see text] and [Formula: see text] that are consistent across multiple simulations, with an average log10-scale standard deviation of 0.28 for on-rates and 0.56 for off-rates, which is well within an order of magnitude and far better than previously observed for previous applications of the WExplore algorithm. Our rank ordering of the three host-guest pairs agrees with the reference calculations, however our predicted [Formula: see text] values were systematically lower than the reference by an average of 4.2 kcal/mol. Using tree network visualizations of the trajectories in the REVO algorithm, and conformation space networks for each system, we analyze the results of our sampling, and hypothesize sources of discrepancy between our [Formula: see text] values and the reference. We also motivate the direct inclusion of [Formula: see text] and [Formula: see text] challenges in future iterations of SAMPL, to further develop the field of ligand binding kinetics prediction and modeling.

Keywords: Binding affinity; Kinetics; Molecular dynamics; REVO; Rare events; SAMPL6; WExplore; Weighted ensemble.

Fig. 1

Structure of the ligands used in this study. Top: quinine, referred to herein as CB8-G3. Middle: 5-hexenoic acid (deprotonated form), referred to herein as OA-G3. Bottom: 4-methyl pentanoic acid (deprotonated form), referred to herein as OA-G6

Fig. 2

The REVO algorithm. Each cycle begins by running an ensemble of walkers forward in time using unbiased dynamics. The distances between the walkers are used to calculate a variance (Eq. 2). In the resampling loop (blue), coupled cloning and merging operations are proposed, and they are accepted only if they result in a higher V. If the proposed V is lower, the resampling loop is terminated and dynamics are continued for the next cycle

Fig. 3

Ensemble splitting. An equilibrium host–guest binding system is split into two nonequilibrium ensembles for the calculation of on and off rates. This is done by defining “bound” and “unbound” basins (left and right of each ensemble). The “unbinding” ensemble (top) is the set of trajectories that have most recently visited the bound basin. The “binding” ensemble (bottom) is the set of trajectories that most recently visited the unbound basin. The on and off rates are directly computed using the time averaged trajectory flux (${\overline{\varphi }}_b$ or ${\overline{\varphi }}_u$) between the ensembles

Fig. 4

Weights of warped walkers. Weights of warping events for the unbinding (top row) and rebinding (bottom row) simulations. In both cases the points are colored according to the index of the corresponding starting pose (0, blue; 1, red; 2, yellow; 3, green; 4, brown)

Fig. 5

Spatial distribution of warped walkers. Structures of warping events for the unbinding simulations viewed from the front and back. Guest ligands are colored according to the index of the corresponding starting pose (0, blue; 1, red; 2, yellow; 3, green; 4, brown)

Fig. 6

Predicted kinetics and free energies. The calculated free energies (top), off-rates (middle), and on-rates (bottom) are shown as a function of simulation time for each starting pose in each host–guest system. The curves are colored according to the index of the starting pose as in Figs. 4 and 5. The calculated binding free energies are compared with experimental measurements (horizontal red line) [38], and the computational reference (dashed black line) for each system

Fig. 7

Trajectory trees show all cloning and merging events in a simulation. The trajectory tree for the first 1329 cycles of the OA-G3–0 unbinding simulation is shown. Each horizontal row in this tree represents a cycle, and the placement of all 48 nodes in the row is determined by minimizing an energy function (see “Visualization of trajectory trees” in “Methods” section). SASA is used to color the nodes, with blue and dark green indicating bound structures, and yellow to orange indicating unbound

Fig. 8

Conformation space networks for the unbinding simulations. Each node in a CSN represents a cluster of host–guest structures. Edges in the networks connect clusters that are seen to interconvert in the REVO simulations. The size of each node is proportional to the number of times it was observed in the unbinding simulations. Nodes are colored according to the solvent accessible surface area of the guest molecule, as shown in the color-bars on the right. The clusters corresponding to the starting poses are labeled in each network