The interaction between a ligand and a protein involves a multitude of conformational states. To achieve a particular deeply-bound pose the ligand must search across a rough free energy landscape, with many metastable minima. Creating maps of the ligand binding landscape is a great challenge, as binding and release events typically occur on timescales that are beyond the reach of molecular simulation. The WExplore enhanced sampling method is well-suited to build these maps, as it is designed to broadly explore free-energy landscapes, and is capable of simulating ligand release pathways that occur on timescales as long as minutes. WExplore also uses only unbiased trajectory segments, allowing for the construction of Markov state models (MSM) and conformation space networks that combine the results of multiple simulations. Here we use WExplore to study two bromodomain-inhibitor systems using multiple docked starting poses (Brd4-MS436 and Baz2B-ICR7), and synthesize our results using a series of MSMs using time-lagged independent component analysis. Ranking the starting poses by exit rate agrees with the crystal structure pose in both cases. We also predict the most stable pose using the equilibrium populations from the MSM, but find that the prediction is not robust as a function of MSM parameters. The simulated trajectories are synthesized into network models that visualize the entire binding landscape for each system, and we examine transition paths between deeply-bound stable states. We find that, on average, transitions between deeply bound states convert through the unbound state 81% of the time, implying a trial-and-error approach to ligand binding. We conclude with a discussion of the implications of this result for both kinetics-based drug discovery and virtual screening pipelines that incorporate molecular dynamics.

[Full text]