How Robust Is the Ligand Binding Transition State?

For many drug targets, it has been shown that the kinetics of drug binding (e.g., on rate and off rate) is more predictive of drug efficacy than thermodynamic quantities alone. This motivates the development of predictive computational models that can be used to optimize compounds on the basis of their kinetics. The structural details underpinning these computational models are found not only in the bound state but also in the short-lived ligand binding transition states. Although transition states cannot be directly observed experimentally due to their extremely short lifetimes, recent successes have demonstrated that modeling the ligand binding transition state is possible with the help of enhanced sampling molecular dynamics methods. Previously, we generated unbinding paths for an inhibitor of soluble epoxide hydrolase (sEH) with a residence time of 11 min. Here, we computationally modeled unbinding events with the weighted ensemble method REVO (resampling of ensembles by variation optimization) for five additional inhibitors of sEH with residence times ranging from 14.25 to 31.75 min, with average prediction accuracy within an order of magnitude. The unbinding ensembles are analyzed in detail, focusing on features of the ligand binding transition state ensembles (TSEs). We find that ligands with similar bound poses can show significant differences in their ligand binding TSEs, in terms of their spatial distribution and protein-ligand interactions. However, we also find similarities across the TSEs when examining more general features such as ligand degrees of freedom. Together these findings show significant challenges for rational, kinetics-based drug design.

Predicting the structural basis of targeted protein degradation by integrating molecular dynamics simulations with structural mass spectrometry.

Targeted protein degradation (TPD) is a promising approach in drug discovery for degrading proteins implicated in diseases. A key step in this process is the formation of a ternary complex where a heterobifunctional molecule induces proximity of an E3 ligase to a protein of interest (POI), thus facilitating ubiquitin transfer to the POI. In this work, we characterize 3 steps in the TPD process. (1) We simulate the ternary complex formation of SMARCA2 bromodomain and VHL E3 ligase by combining hydrogen-deuterium exchange mass spectrometry with weighted ensemble molecular dynamics (MD). (2) We characterize the conformational heterogeneity of the ternary complex using Hamiltonian replica exchange simulations and small-angle X-ray scattering. (3) We assess the ubiquitination of the POI in the context of the full Cullin-RING Ligase, confirming experimental ubiquitinomics results. Differences in degradation efficiency can be explained by the proximity of lysine residues on the POI relative to ubiquitin.

Creating Maps of the Ligand Binding Landscape for Kinetics-Based Drug Discovery.

Simulations of ligand-protein interactions can be very useful for drug design and to gain biological insight. Full pathways of ligand-protein binding can be used to get information about ligand binding transition states, which form the rate-limiting step of the binding and release processes. However, these simulations are typically limited by the presence of large energy barriers that separate stable poses of interest. Here we describe a simulation protocol for exploring and analyzing landscapes of ligand-protein interactions that makes use of molecular docking, enhanced molecular simulation with the weighted ensemble algorithm, and network analysis. It can be accomplished using a modest cluster of graphics processing units and freely accessible software. This protocol focuses on the construction and analysis of a network model of ligand binding poses and provides links to resources that describe the other steps in more detail. The end result of this protocol is a map of the ligand-protein binding landscape that identifies transition states of the ligand binding pathway, as well as alternative bound poses that could be stabilized with modifications to the ligand.

Wepy: A Flexible Software Framework for Simulating Rare Events with Weighted Ensemble Resampling.

Here, we introduce the open-source software framework wepy (https://github.com/ADicksonLab/wepy) which is a toolkit for running and analyzing weighted ensemble (WE) simulations. The wepy toolkit is in pure Python and as such is highly portable and extensible, making it an excellent platform to develop and use new WE resampling algorithms such as WExplore, REVO, and others while leveraging the entire Python ecosystem. In addition, wepy simplifies WE-specific analyses by defining out-of-core tree-like data structures using the cross-platform HDF5 file format. In this paper, we discuss the motivations and challenges for simulating rare events in biomolecular systems. As has previously been shown, high-dimensional WE resampling algorithms such as WExplore and REVO have been successful at these tasks, especially for rare events that are difficult to describe by one or two collective variables. We explain in detail how wepy facilitates implementation of these algorithms, as well as aids in analyzing the unique structure of WE simulation results. To explain how wepy and WE work in general, we describe the mathematical formalism of WE, an overview of the architecture of wepy, and provide code examples of how to construct, run, and analyze simulation results for a protein-ligand system (T4 Lysozyme in an implicit solvent). This paper is written with a variety of readers in mind, including (1) those curious about how to leverage WE rare-event simulations for their domain, (2) current WE users who want to begin using new high-dimensional resamplers such as WExplore and REVO, and (3) expert users who would like to prototype or implement their own algorithms that can be easily adopted by others.

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

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.

Unbiased Molecular Dynamics of 11 min Timescale Drug Unbinding Reveals Transition State Stabilizing Interactions.

Ligand (un)binding kinetics is being recognized as a determinant of drug specificity and efficacy in an increasing number of systems. However, the calculation of kinetics and the simulation of drug unbinding is more difficult than computing thermodynamic quantities, such as binding free energies. Here we present the first full simulations of an unbinding process at pharmacologically relevant time scales (11 min), without the use of biasing forces, detailed prior knowledge, or specialized processors using the weighted ensemble based algorithm, WExplore. These simulations show the inhibitor TPPU unbinding from its enzyme target soluble epoxide hydrolase, which is a clinically relevant target that has attracted interest in kinetics optimization in order to increase efficacy. We make use of conformation space networks that allow us to conceptualize unbinding not just as a linear process, but as a network of interconnected states that connect the bound and unbound states. This allows us to visualize patterns in hydrogen-bonding, solvation, and nonequilibrium free energies, without projection onto progress coordinates. The topology and layout of the network reveal multiple unbinding pathways, and other rare events, such as the reversal of ligand orientation within the binding site. Furthermore, we make a prediction of the transition state ensemble, using transition path theory, and identify protein-ligand interactions which are stabilizing to the transition state. Additionally, we uncover trends in ligand and binding site solvation that corroborate experimental evidence from more classical structure kinetics relationships and generate new questions as to the role of drug modifications in kinetics optimization. Finally, from only 6 µs of simulation time we observed 75 unbinding events from which we calculate a residence time of 42 s, and a standard error range of 23 to 280 s. This nearly encompasses the experimental residence time 11 min (660 s). In addition to the insights to sEH inhibitor unbinding, this study shows that simulations of complex processes on timescales as long as minutes are becoming feasible for more researchers to perform.

Genetic basis of photosynthetic responses to cold in two locally adapted populations of Arabidopsis thaliana.

Local adaptation is common, but the traits and genes involved are often unknown. Physiological responses to cold probably contribute to local adaptation in wide-ranging species, but the genetic basis underlying natural variation in these traits has rarely been studied. Using a recombinant inbred (495 lines) mapping population from locally adapted populations of Arabidopsis thaliana from Sweden and Italy, we grew plants at low temperature and mapped quantitative trait loci (QTLs) for traits related to photosynthesis: maximal quantum efficiency (Fv/Fm), rapidly reversible photoprotection (NPQfast), and photoinhibition of PSII (NPQslow) using high-throughput, whole-plant measures of chlorophyll fluorescence. In response to cold, the Swedish line had greater values for all traits, and for every trait, large effect QTLs contributed to parental differences. We found one major QTL affecting all traits, as well as unique major QTLs for each trait. Six trait QTLs overlapped with previously published locally adaptive QTLs based on fitness measured in the native environments over 3 years. Our results demonstrate that photosynthetic responses to cold can vary dramatically within a species, and may predominantly be caused by a few QTLs of large effect. Some photosynthesis traits and QTLs probably contribute to local adaptation in this system.

Multiple Ligand Unbinding Pathways and Ligand-Induced Destabilization Revealed by WExplore.

We report simulations of full ligand exit pathways for the trypsin-benzamidine system, generated using the sampling technique WExplore. WExplore is able to observe millisecond-scale unbinding events using many nanosecond-scale trajectories that are run without introducing biasing forces. The algorithm generates rare events by dividing the coordinate space into regions, on-the-fly, and balancing computational effort between regions through cloning and merging steps, as in the weighted ensemble method. The averaged exit flux yields a ligand exit rate of 180 µs, which is within an order of magnitude of the experimental value. We obtain broad sampling of ligand exit pathways, and visualize our findings using conformation space networks. The analysis shows three distinct exit channels, two of which are formed through large, rare motions of the loop regions in trypsin. This broad set of ligand-bound poses is then used to investigate general properties of ligand binding: we observe both a direct stabilizing effect of ligand-protein interactions and an indirect destabilizing effect on intraprotein interactions that is induced by the ligand. Significantly, the crystallographic binding poses are distinguished not only because their ligands induce large stabilizing effects, but also because they induce relatively low indirect destabilizations.

Ligand Release Pathways Obtained with WExplore: Residence Times and Mechanisms.

The binding of ligands with their molecular receptors is of tremendous importance in biology. Although much emphasis has been placed on characterizing binding sites and bound poses that determine the binding thermodynamics, the pathway by which a ligand binds importantly determines the binding kinetics. The computational study of entire unbiased ligand binding and release pathways is still an emerging field, made possible only recently by advances in computational hardware and sampling methodologies. We have developed one such method (WExplore) that is based on a weighted ensemble of trajectories, which we apply to ligand release for the first time, using a set of three previously characterized interactions between low-affinity ligands and the protein FKBP-12 (FK-506 binding protein). WExplore is found to be more efficient that conventional sampling, even for the nanosecond-scale unbinding events observed here. From a nonequilibrium ensemble of unbinding trajectories, we obtain ligand residence times and release pathways without using biasing forces or a Markovian assumption of transitions between regions. We introduce a set of analysis tools for unbinding transition pathways, including using von Mises-Fisher distributions to model clouds of ligand exit points, which provide a quantitative proxy for ligand surface diffusion. Differences between the transition pathway ensembles of the three ligands are identified and discussed.