Ligand Release Pathways Obtained with WExplore: Residence Times

Subscriber access provided by UCL Library Services

Article

Ligand Release Pathways Obtained with WExplore: Residence Times and Mechanisms Alex Dickson, and Samuel D Lotz J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b04012 • Publication Date (Web): 27 May 2016 Downloaded from http://pubs.acs.org on May 31, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Ligand Release Pathways Obtained with WExplore: Residence Times and Mechanisms Alex Dickson∗,†,‡ and Samuel D. Lotz† Department of Biochemistry & Molecular Biology, Michigan State University, East Lansing, MI, and Department of Computational Mathematics, Science and Engineering, Michigan State University, East Lansing, MI E-mail: [email protected]\\Phone:+1517-884-8985

∗

To whom correspondence should be addressed Biochemistry ‡ CMSE †

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract 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, that provide a quantitative proxy for ligand surface diffusion. Differences between the transition pathway ensembles of the three ligands are identified and discussed.

Introduction Computational studies of ligand binding are dominated by methods that consider only the endpoints of the process. Docking is widely used to determine bound ligand poses and to estimate binding free energies using empirical scoring functions, and can be conducted on a large scale to facilitate receptor-based virtual screening. 1–3 More robust methods to calculate binding free energies also typically consider only the bound and unbound states of the ligand. 4–7 Accordingly, while much is known about the relationship between ligand structure

2


Page 2 of 32

Page 3 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and binding affinity, much less is known about the general properties of the pathways of small molecule binding. Simulations employing “steering” forces to encourage ligands to unbind from their pockets 8–11 are useful to determine general properties of pathways, especially of deeply buried ligands, but they can introduce significant perturbations to unbinding mechanisms. Other techniques such as metadynamics 12–14 and umbrella sampling 15,16 can recapitulate binding pathways as well as binding kinetics, 17,18 although a good choice of collective variables is crucial, and often not obvious. Enabled by advances in computer hardware and sampling algorithms, full unbiased ligand binding trajectories have emerged for a number of systems. This includes long, continuous trajectories of ligand binding, 19 as well as large numbers of shorter trajectories synthesized by Markov state models. 20–22 It was shown recently that a dramatic sampling improvement could be obtained by iteratively constructing Markov state models and using these to inform sampling on-the-fly, with trajectory management steps happening on the 10 ns timescale. 23 We have developed a similar method, WExplore, 24 which also uses unbiased trajectories that are actively managed on a much shorter timescale (20 ps), but does not rely on a Markovian assumption. WExplore defines a set of regions in a high-dimensional order parameter space, and uses trajectory cloning and merging steps – as in the weighted ensemble algorithm 25 – to improve the sampling of low-probability regions, such as transition states that lie between two high-probability basins. In this work we apply WExplore to protein-ligand systems for the first time, allowing for a more thorough exploration of ligand unbinding transition pathways, which, according to the time-reversibility principle, reveal the binding transition paths as well. Previous ligand binding simulations have provided a mixed picture of general ligand binding mechanisms, particularly in regard to the extent of surface diffusion prior to ligand release. Trajectories of dasatinib binding to Src kinase from Shan et al 19 have shown extensive surface diffusion prior to binding, involving almost the entire surface of the protein (this behavior is referred to herein as the “surface diffusion model”). Moderate surface diffusion

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

was found trajectories of benzamidine binding to trypsin. 20 In contrast, some enhanced sampling restraint potentials, such as funnel metadynamics, 14 do not allow for substantial surface diffusion in their predictions of the ligand binding free energy maps. As these reconstructed binding pathways establish first contact with the protein directly at the binding site, we refer to this as the “dartboard model”. These two models offer very different pictures of ligand binding processes, and to distinguish between them we need to develop general quantitative metrics of surface diffusion on ligand binding transition pathways. In this paper we define such a quantitative measure by collecting the positions of the ligand as it exits the environment of the protein. These exit point clouds can be analyzed using probability maps that are projected onto a unit sphere, the center of which is mapped to the center of mass of the protein. The concentration of points on this sphere can be modeled quantitatively using Von Mises-Fisher probability distributions, 26,27 which are characterized by a mean vector and a concentration parameter κ. This provides a direct proxy of surface diffusion: for systems that exhibit high surface diffusion during binding and unbinding pathways, the exit points will be very disperse, with low κ values; for systems with minimal diffusion, the κ values will be very high. This general strategy allows for robust comparison both between ligands binding the same protein, and between different ligand-protein pairs. We use the WExplore method to examine ligand release distributions for a well-studied ligand binding model system with three ligands: FK506 binding protein (referred to here as FKBP) with three low affinity ligands. FKBP belongs to the peptidylprolyl cis/trans isomerase (PPIases) family of proteins. 28,29 It has been the subject of much interest since its discovery, with structures and free energies of binding determined to a large number of drug-like ligands, both experimentally 29–34 and computationally. 35 The interaction of FKBP with a set of low affinity small-molecule ligands has also been characterized, with bound structures and biophysical properties. 36 Further analysis was conducted recently in silico by Huang and Caflisch, 21 which determined the landscape of binding and unbinding pathways in great detail, and determined the unbinding kinetics for a series of ligands. This extensive

4


Page 4 of 32

Page 5 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


characterization, as well as the fast unbinding time that is tractable to conventional molecular simulation, makes FKBP an ideal test case to demonstrate WExplore sampling applied to ligand binding systems for the first time.

Methods Molecular dynamics sampling The trajectory segments are run using CHARMM 37 through the OpenMM-CHARMM interface, with the CHARMM36 force field for the protein. Ligand parameterization is done using the CHARMM Generalized Force Field (CGenFF). 38,39 Solvation, minimization and equilibration are done separately for the three ligands BUT, DMSO and DSS, using structures 1D7J, 1D7H and 1D7I, respectively 36 (Figure 1). Each structure is solvated in a cubic box allowing 12 ˚ A of space between the protein and the cell boundaries, resulting in 13623, 13512 and 12350 water molecules for BUT, DMSO and DSS, respectively. One chloride ion is used to neutralize each system. We use SHAKE constrain the hydrogen atoms along their covalent bonds, with a tolerance of 10−8 . Particle-mesh Ewald summation is used for nonbonded interactions, along with a van der Waals switching function that scales the nonbonded interactions to zero at 10 ˚ A, starting from 8.5 ˚ A.

BUT

DMSO

DSS

Figure 1: Three low-affinity ligands used in this study. Left: 4-hydroxy-2-butanone is referred to throughout as “BUT”. Middle: Dimethylsulfoxide is referred to as “DMSO”. Right: Methyl sulfinyl-methyl sulfoxide is referred to as “DSS”. The system is minimized with constraints on solute atoms for 500 steps using steepest 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

descent and 500 steps using the adopted basis Newton Raphson method. The constraints are removed and the minimization is repeated. The system is then heated from 50 K to 300 K over 100 ps, using 2 fs timesteps, followed by equilibration for 0.5 ns at 300 K. The resulting structure is used as the initial coordinates for all replicas in the WExplore method. Langevin dynamics is performed at 300 K with a friction coefficient of 1.0 ps−1 . We use a constant pressure ensemble with a Monte Carlo barostat coupled to a reference pressure of 1 atm, and volume moves attempted every 50 steps. We use a group of eight graphics processing units, each running their own instance of CHARMM with OpenMM 6.3 40 to run trajectories for the replicas. These trajectories are managed by a script that implements the WExplore method by communicating with the CHARMM instances on-the-fly.

WExplore sampling The WExplore method 24 is an extension of the weighted ensemble (WE) algorithm 25 that facilitates usage in high-dimensional spaces, and is particularly useful in simulating lowentropy to high-entropy transitions, such as protein unfolding or ligand release. In the original WE algorithm, sampling of low-probability states is enhanced by defining a set of regions that span the space, and enforcing that all regions are sampled evenly. Multiple copies of the system (“walkers”) are run in parallel (here, 48), and each is given a weight with which it contributes to statistical averages; the sum of weights over all walkers is equal to 1. Periodically throughout the simulation, walkers are merged in over-represented regions and cloned in under-represented regions in order to keep the number of walkers in each region as even as possible. When a walker is cloned, its weight is split up evenly amongst the clones. When two walkers A and B, with weights wA and wB , are merged, the resulting walker has weight wA + wB , and takes on the configuration of walker A with probability wA /(wA + wB ), and otherwise takes on the configuration of walker B. WExplore is an extension to WE that allows for efficient application to high-dimensional spaces, and allows for the number of walkers to be far less than the number of regions. A 6


Page 6 of 32

Page 7 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


more thorough description of the algorithm can be found in previous work. 24,41 Here we discuss the main details of the algorithm, focusing on aspects that are unique to this study. WExplore uses Voronoi polyhedra to define sampling regions in a high-dimensional space. These Voronoi polyhedra are defined by a set of “images”, which are conformations of the system, and are defined to all be a critical distance apart (dc ) from each other, using an RMSD-based (root mean squared distance) distance metric. Importantly, these regions are defined using a hierarchy: the top level regions tile the entire space, and these are tiled by smaller regions, which are themselves tiled by even smaller regions, and so on. Each level of the hierarchy has its own value of dc , with the top level regions having the highest, and the bottom having the lowest. This hierarchical structure allows for replica balancing across multiple length scales: a balancing procedure is first performed at the top level of the hierarchy, and then at each lower level, prioritizing balancing computational effort amongst regions that are the farthest apart. The structure also helps with region assignment: at each level of the region “tree”, the distance to each image needs to be calculated for only a small number of branches. Here we apply WExplore to study small molecule binding and release for the first time. As noted above, the key quantity to run a WExplore simulation is a way of measuring the distance between two conformations of the system. Our distance metric for small molecule binding is defined as follows. First the two conformations are aligned using a selection of residues that are local to the binding site. These are determined as the set of protein residues that are within a cutoff of 8 ˚ A in the crystal structure binding pose. The distance is then determined using the RMSD of the absolute positions of the ligand, that is without any further translation or rotation. This metric captures both internal ligand motions, as well as ligand translations and rotations. In WExplore, regions are defined dynamically over the course of the simulation, and do not require any a priori information about the dynamics of interest. We use a 4-level sampling hierarchy with region sizes of 2.5, 3.5, 5.0 and 10 ˚ A, and a maximum number of

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 32

images of 10 at each sampling resolution, for a theoretical maximum of 10000 regions. In a sampling cycle, all 48 walkers are run forward in time for 20 ps, each is assigned to a region, and trajectories are cloned and merged according to the balancing algorithm. 24 The aggregate simulation time of a WExplore simulation after 250 cycles is then 240 ns. We run two WExplore simulations for each of the three ligands, as well as a conventional simulation of the same length (48 trajectories, 240 ns total). In total 2.16 µs of simulation is performed. Two aspects of our implementation are significantly different in previous work. Firstly, to minimize the sampling expenditure on conformations where the ligand is far away from the protein, we apply the following rule: if a walker has a minimum interatomic distance to the protein that is between 5 and 10 ˚ A, it is not allowed to define new regions, and it is excluded from walker cloning and merging operations. Secondly, we institute a minimum and a maximum weight for the walkers, and we do not clone walkers that are below the minimum weight (pmin ), and do not merge walkers that are above the maximum weight (pmax ). pmax is set to be 0.1, in order to ensure that there is always at least 5 to 10 walkers in the original binding site. In order to efficiently determine unbinding rates, pmin should ideally be set close to the probability of the rarest conformations that are still of importance to the unbinding ensemble. In other words, pmin should be close to the probability of conformations in the transition state, but since this is not known a priori, we use the following procedure. We initialize pmin at pimin = 10−40 , and allow this to grow as follows:

pmin =

    pimin    

if poff /10 < pimin

poff /10 if pimin < poff /10 < 1/2N        1 otherwise 2N

(1)

where poff is the sum of the weights of all exited trajectories collected so far and N is the number of walkers in the simulation. The net effect of this procedure is that even after we have an estimate of the transition state probability, we are still giving the system the

8


Page 9 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


resources it needs to make transitions with higher transition state probabilities.

Exit rates Following previous work using Nonequilibrium Umbrella Sampling 42,43 and the string method, 44 we can determine the rate of ligand release in a WExplore simulation by examining the probability flux into a given basin of attraction, averaged over a period of time. For our purposes, consider basin A as the set of ligand configurations that are within a certain RMSD cutoff of the native binding pose (e.g. 3 ˚ A). Consider basin B as the set of ligand configurations where the closest protein-ligand interatomic distance is greater than another cutoff, dU , set here to 10 ˚ A. Here, we simulate the nonequilibrium ensemble where trajectories are initialized in basin A and are terminated in basin B, and the exit rates are formally defined as the transition rate from state A to B. In practice, every cycle we check whether any trajectories have crossed into basin B. For those that have, we record those conformations for analysis, add the weight of that trajectory to a running total, and allow the conformation of that walker to be overwritten by another walker in the next cloning and merging step. At any time, the rate can be calculated by dividing the running total by the total simulation time elapsed up to that point.

Von Mises-Fisher exit distributions Our major goal of this work is to model the collection of ligand exit points and quantitatively determine their variance, or spread. We take advantage of existing tools from directional statistics to accomplish this goal. Von Mises-Fisher (VMF) distributions model probability distributions of points arranged on spheres that are characterized only by a mean direction µ and a concentration parameter κ. In three dimensions, the probability density function is given by f (X) =

κ T eκµ X 2π sinh(κ)

9


(2)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32

where X is any vector on the unit sphere. The exponential term penalizes vectors that are not close to the average vector µ, and κ scales the strength of that penalization. Figure 2 shows representative VMF distributions with 300 points each that are generated using the modified Ulrich’s algorithm. 45 When the concentration parameter κ is large, the spread of points on the sphere is small. For values as small as κ = 1 the distribution of points is almost isotropic on the sphere. front

back

20

>4

5

3 2 1

1

0 -ln(P/Pmax)

Figure 2: Von Mises-Fisher distributions. Random vector sets are shown that correspond to Von Mises-Fisher distributions with varying values of the concentration parameter κ. The spheres are oriented with the average vector aligned along the y-axis, which points out of the page. The patches on the sphere are colored according to the negative logarithm of the probability, with white showing the most probable regions and black the least probable. To visualize our distributions of exit points, we need to first define a common alignment frame, to facilitate comparison of distributions from one simulation to another. We first translate the system to place the center of mass of the protein at the origin. We then rotate the system to orient the center of mass of residues in the binding pocket (defined as those residues that are within 8 ˚ A of the ligand in its crystal structure) along the y-axis. Another 10


Page 11 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


“capping” residue of the protein is arbitrarily chosen (here, FKBP residue 181), and the system is rotated around the y-axis until the projection of the center of mass of the capping residue onto the xz-plane is aligned to the z-axis. After this alignment procedure, we visualize the distribution using a histogram with bins defined using spherical polar coordinates φ and ψ. The plots used here use 30 bins along each coordinate. For the WExplore simulations, each ligand exit point contributes to a histogram according to the weight of the trajectory that recorded the point. For each exit point distribution we obtain the concentration parameter κ using the approximation: 27 κ= where R =

||

P

xi || , n

3R − R3 1 − R2

(3)

and xi is the normalized ith exit point vector, and n is the number

of exit points recorded. We have validated this approximation using our generated VMF distributions above, and found it is accurate for our range of parameters.

Clustering Clustering is performed using MSMBuilder version 3.3.1, 46 using the “RegularSpatial” algorithm. 47 For each ligand we use a set of protein-ligand distances connecting three points in the ligand to the 50 closest heavy atoms in the protein: 150 distances total. We use the Canberra distance metric between the two sets of distances:

d(A, B) =

X |ai − bi | |ai | + |bi | i

(4)

where a is the set of atomic distances in structure A, and b is the set of atomic distances in structure B. The Canberra metric is appropriate for our purposes as it emphasizes differences in distances that are small (i.e., close to the binding pocket). A cluster size of 10 is used in each case, resulting in 485, 497 and 565 cluster centers for BUT, DMSO and DSS, respectively. 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 32

Results & Discussion Residence times As this is the first application of WExplore to a protein-ligand system, we used conventional simulations as a benchmark, both those conducted here and simulations run previously. 21 The final estimates of the mean first passage time (MFPT) of unbinding are shown in Table 1, and show solid quantitative and qualitative agreement between WExplore and conventional simulations conducted here. Averaging over the three simulations: BUT shows the fastest dissociation at 5 ns, DMSO is the next fastest at 8 ns, and DSS is the slowest with a MFPT of 23 ns. Previous simulations documented similar results, 21 and the discrepancies are comparable to estimates of standard error. Interestingly, experimental binding affinities only partially correlate with these results. 36 DSS has the highest binding affinity, but only by a factor of 2. The binding affinity of DMSO was found to be 400 times weaker than BUT (and 800 times weaker than DSS), however the residence times of DMSO and BUT are comparable. Table 1: Mean first passage times of unbinding for three ligands from FKBP, obtained by WExplore simulation (two replicates), straight-forward simulation, and previous simulations by Huang et. al. 21 In the first three cases, standard errors are obtained by block averaging the last 50% of the data into 10 evenly spaced blocks, in each case. The experimentally-determined binding affinity is also shown for each molecule.

BUT DMSO DSS

WExplore1 4±1 5±1 19 ± 5

Unbinding time [ns] WExplore2 SF Huang 2011 7±2 4±1 8±2 13 ± 3 6±1 4±1 27 ± 6 22 ± 4 18 ± 3

Binding affinity [µM] Experiment 1 500 20000 250

Figure S1 shows the convergence of the MFPT estimates as a function of time for both WExplore simulations and the conventional simulations. The estimates from WExplore simulations converge to their final value after approximately 100 ns of total simulation, indicating that 240 ns is adequate to determine unbinding kinetics in these systems. Although 12


Page 13 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


all three ligands dissociate quickly enough to be captured by conventional simulation, the slowest dissociating ligand (DSS) clearly shows the benefits of the WExplore method. For conventional simulation the first exit point is not observed until after over 104 ns of total simulation, while exit paths are observed in WExplore simulations using 13.4 and 17.3 ns of total simulation, respectively. This is encouraging for the future performance of WExplore simulations for high-affinity ligands, since the efficiency gains of weighted ensemble sampling have been shown to increase with the height of the free energy barrier. 48–51

Exit point distributions Even for low affinity ligands, WExplore simulation generates distributions of exit points much more efficiently than conventional MD. Figure 3 shows the number of exit points generated as a function of simulation time for both WExplore (the average of two simulations) and conventional MD. In each case, the most points are determined for BUT, the second most for DMSO, and the least points were determined for DSS. This ordering is what would be expected from the exit rates for the three ligands. For DSS this is particularly stark: only 9 exit points are obtained in 250 ns of simulation with conventional sampling, compared to 101 and 95 from the two WExplore simulations. Figure 4 shows the complete set of exit points determined in the WExplore and conventional sampling simulations for the three ligands DMSO, BUT and DSS. These are all aligned using a common reference frame to facilitate a direct comparison between maps. All maps show a high degree of variation, indicating substantial diffusion both at and near the surface of the protein. It is important to note that these points represent the first time in that trajectory that the ligand has been more than 10 ˚ A away from the surface of the protein, as measured by the closest atom-atom distance; throughout the unbinding trajectories the ligands are closer than 10 ˚ A. We quantify this spread by fitting these exit point clouds to von Mises-Fisher distributions, which are characterized by both a mean vector and a concentration parameter, κ. The concentration parameters found here for these three ligands 13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Number of exit points


200 180 160

Page 14 of 32

BUT

WExplore Conventional

140 120 100

DMSO DSS

80 60 40 20 0 0

50

100

150

200

250

Simulation time (ns)

Figure 3: WExplore more efficiently generates exit points. A running total of the number of exit points observed over the course of the simulation is plotted for both WExplore (thick lines) and conventional sampling (thin lines). The WExplore curves are the averages of the two WExplore simulations performed for each ligand. are relatively close: BUT is only slightly less concentrated than the other two ligands. The absolute values of the κ parameters are also telling. The values here (κ ≈ 5) are intermediate between what would be expected in a “surface crawling” model and a “dartboard” model (as described in the Introduction). Convergence of the measured κ values for the individual simulations are shown in Figure S2. We observe that different WExplore simulations can converge to slightly different values in κ, which suggests better results can come from averaging over a set of multiple WExplore simulations. The mean vectors of escape, as used in the VMF distributions for each ligand, are shown oriented onto the protein structure in Figure 5. This clearly reveals differences in the BUT ensemble compared to the other two ligands. The mean vector of escape is rotated toward the lower left-hand quadrant, which is also shown in the probability maps for both WExplore and conventional sampling in Figure 4. As shown in the side view of Figure 5, the protein is not smooth along its front face, but has a bulge in the lower left. As all three mean escape vectors point away from the bulge, it is likely that this acts as a steric barrier that drives

14


Page 15 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A

C

WExplore

Conventional sampling 7

DMSO κ=5.8±0.9 0 front

back

front

back 7

B

BUT κ=4.0±0.2 0 front

back

front

back 7

DSS κ=5.5±0.9 front

back

front

0 back -ln(P/P ) max

Figure 4: Mapping ligand exit distributions. (A) The FKBP protein is shown, oriented as described in Methods: the vector connecting the center of mass of the protein and the center of mass of the binding pocket residues is aligned to the y-axis, which points out of the page. This orientation can be used as a reference to understand panels (B) and (C). The equilibrated bound pose of the ligand DMSO is shown in CPK representation. (B) The cloud of exit points from a WExplore simulation of DMSO. Note that each of these conformations contributes to the spherical maps in (C) with a different statistical weight. (C) Spherical heat maps of exit point distributions for WExplore and conventional sampling. For each ligand, the front and back of the sphere are shown side-by-side. The front view corresponds to the view shown in (A), and the back corresponds to a 180◦ rotation around the vertical axis of that structure. Each spherical map shows the free energy landscape of exit points projected onto the unit sphere. The most probable exit location is shown in white, and the least probable is shown in a dark color (DMSO = green, BUT = red, DSS = purple). The κ values shown are the averages of the two WExplore simulations and the conventional simulations, the uncertainty is the standard error across all three simulations.

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

exit pathways to the right.

DMSO DSS BUT

90°

Figure 5: Mean escape vectors. Mean vectors of the von Mises-Fisher probability distributions are displayed on the FKBP protein. In the left frame the protein is oriented as in Figure 4, and the right frame shows a side view. DMSO and DSS are both pointing slightly up and to the right, while BUT is rotated significantly toward the bottom left quadrant.

Surface binding maps To help investigate the origins of the difference between BUT and the other two ligands we examined ligand binding to various regions of the protein along the exit pathways observed in the WExplore simulations for each ligand. Figure 6 compares binding probabilities, as measured using an atom-atom cutoff distance of 6 ˚ A, for the three different ligands mapped onto the protein surface. The views on the left hand side show the front of the molecule, with the binding pocket in the middle. In accord with the results above, BUT shows a much higher density of bound ligand in the lower left hand quadrant than either DSS or DMSO. This can be seen more clearly in the side views for each ligand, shown in the right hand column. The surface density for DSS is much patchier than for the other two ligands, indicating a smaller number of more tightly bound poses. To investigate the structural details of these poses we perform a clustering procedure and pull out representative conformers (see Methods). The same cluster radius is used for each ligand allowing for comparison across the ligand set. We define a pseudo Shannon

16


Page 16 of 32

Page 17 of 32


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

entropy as S = −kT

P

i

pi ln pi , where the sum is over the set of clusters, and the probability

of a cluster (pi ) is the sum of the probabilities of each conformation in the cluster. Note that this is not equal to the true Shannon entropy, as we are only using statistics from the nonequilibrium unbinding ensemble. Table 2 shows the S values for the three ligands, and these are consistent with the ligand surface densities: DSS has a significantly lower conformational entropy in the bound ensemble. Interestingly, this is not reflected in the number of clusters that are found for DSS, which underscores the ability of WExplore to broadly sample low probability conformational ensembles. Table 2: Pseudo-Shannon entropy of the bound ensemble, and the number of clusters found for each ligand, using a constant cluster radius. DMSO BUT DSS

Nc 497 485 565

S (kT) 4.31 4.49 3.59

Binding pose analysis The most probable bound structures are also determined by choosing a representative structure from the highest probability regions for each ligand. These structures are overlaid with the crystallographic poses in Figure S3. After aligning the structures based on a set of FKBP residues in the binding pocket, the root mean squared distance between the highest probability poses and the crystallographic poses is determined. All three ligands show general agreement with respect to the position of the ligand. For DSS and DMSO, these poses show high agreement with the crystal structure poses, with RMSD values of 1.03 ˚ A and 0.77 ˚ A respectively. The highest probability BUT pose found here differs in orientation from the pose in the crystal structure with PDB ID 1D7J 36 (RMSD = 3.08 ˚ A), with the carbonyl group pointing the opposite direction (in our case, towards the NH group of residue I56). Interestingly, this pose was also identified in a set of high probability poses of BUT identified by Huang and Caflisch, 21 although it was not found to be the most probable. 18


Page 18 of 32

Page 19 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Noncovalent interactions for hydrogen bonding and pi-cation interactions are then profiled with the PLIP (Protein-Ligand Interaction Profiler) program. 52 For hydrogen bonds to be detected the acceptor and donor must be less than 4.1 ˚ A apart and at an angle between 180◦ and 100◦ . 53 Pi-cation interactions are detected when the charged atom and aromatic ring center are less than 6.0 ˚ A apart. 54 We identify interactions in the representative conformers of high-weight cluster centers (w > 0.005), and show all such unique interactions for each ligand in Figure 7. BUT forms hydrogen bonding interactions in 128 representative structures (compared to 88 and 69 for DSS and DMSO, respectively), contributing to its higher entropy of clustering and its broad surface density. In contrast, DSS has specific high-probability pi-cation interactions between the partially positive sulfoxide (0.310 D) and the aromatic rings of TRP59, specifically the −0.61 D indole N (N1) (panel DSS1 of Figure 7) and the −0.115 D carbons of the six member indole ring (C4-7) (panel DSS2 of Figure 7). The interaction with N1 has a much higher probability (w = 0.274) than the most probable interactions for BUT (w = 0.088) and DMSO (w = 0.042). The higher specificity interactions in DSS are consistent with its more restricted surface densities (Figure 6), and lower entropy (Table 2) than both BUT and DMSO. It is interesting that even though DMSO is a structural subset of DSS (also with a 0.31 D sulfoxide), we did not observe probable poses involving pi-cation interactions nor were any poses stabilized to the extent of DSS1. This is likely due to the small size of DMSO, which lacks the steric confinement that can stabilize electrostatic interactions in larger ligands like DSS. This suggests that fragment-based screening approaches using very small ligands may miss important electrostatic interactions such as pi-cation interactions, salt-bridges, and halogen bonds that can help drive receptor specificity.

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

BUT1 w=0.088

BUT2 w=0.017 Asp37

DSS1 w=0.274

DSS2 w=0.037

Asp37

Tyr26 BUT

Tyr82

DMSO1 w=0.042

Page 20 of 32

DSS

DMSO

Trp59

Tyr82

Glu54 Ile56 BUT3 w=0.011

Ile56

Ile56

BUT4 w=0.008

DMSO2 w=0.010

DSS3 w=0.010

DSS4 w=0.008

Figure 7: High probability ligand-protein interactions. Representative structures for cluster centers with unique ligand-protein interactions and weight greater than 0.005 are shown for each ligand numbered in descending order of cluster weight. Red dashed lines indicate hydrogen bonds and solid blue lines indicate pi-cation interactions.

Conclusion The study marks, to our knowledge, the first analysis of ligand exit distributions. The approach presented here for their visualization and analysis is general, quantitative and can facilitate comparison between different ligands and different molecular receptors as well. We thus expect that concentration parameter (κ) values and unit sphere projections as in Figure 4 will be useful proxies to quantify ligand surface diffusion prior to release for a large number of ligands, including larger drug-like molecules. Of particular interest will be the dependence of κ on ligand properties such as binding affinity and ligand hydrophobicity. Also, for larger proteins we expect larger κ values, as it becomes harder for a ligand to diffuse from one side of a protein to the other before diffusing away from the protein surface. Membrane proteins also present a special case, where the membrane provides a barrier to ligand surface diffusion. Ligand release in general is an extremely demanding computational problem, and even for these low affinity ligands, the enhanced sampling method WExplore enabled our collection of a statistically relevant set of exit pathways. This work is the first demonstration of WExplore on a ligand binding system. The distance metric used to define the sampling regions captures 20


Page 21 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


both internal ligand motions as well as global ligand translations and rotations with respect to the binding pocket. We expect that this distance metric will also work well for flexible ligands with binding pathways that feature motion along internal degrees of freedom. One limitation is that this distance metric does not enhance sampling of protein motions that are distant from the binding site. Thus if ligand release is triggered by allosteric motions far from the binding site, a different distance metric should be used that would also incorporate those motions in order to efficiently sample ligand release. WExplore is well suited to ligand binding applications as it does not require a specific single order parameter, and thus does not bias our collection of exit paths in any particular direction. Surface diffusion estimates from WExplore are thus more reliable than those obtained using pulling forces, for example. Numerous other undirected enhanced sampling methods such as temperature accelerated molecular dynamics, 55,56 orthogonal space random walk, 57 self-guided Langevin dynamics, 58 random acceleration molecular dynamics 59 and random expulsion molecular dynamics 60 are also suitable for this purpose, although WExplore differs in that trajectories are run using the unbiased Hamiltonian, allowing for direct incorporation into Markov state modeling frameworks 61 and network visualization tools. 62,63 Looking forward, a robust comparison of the performance of enhanced sampling methods for the generation of ligand binding and unbinding trajectories is needed, from the standpoint of efficiency as well as accuracy. The coming years should see a growth in simulations of full binding and unbinding trajectories, due to growing computational power as well as improvements in ligand force field parameterization. General, quantitative metrics that measure properties of ligand binding pathways, together with visualization tools will go a long way to help define principles and models for the basic physics of ligand binding. A deep physical understanding of the binding process will improve our understanding of binding kinetics, helping both directly to design drugs that have the desired kinetic properties – such as long residence times 64–66 – and indirectly to inform scoring functions used in ligand- and receptor-based virtual screening.

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Acknowledgement A.R.D. would like to thank Drs. Jose Parea and Yuying Xie for helpful discussions. We also acknowledge support from the High Performance Computing Center at Michigan State University.

Supporting Information Available Figure S1: Estimated mean first passage times for unbinding as a function of simulation time. Figure S2: Estimated concentration parameter values as a function of simulation time. Figure S3: Comparison of highest probability poses with crystallographic poses.

References (1) Cheng, T.; Li, Q.; Zhou, Z.; Wang, Y.; Bryant, S. H. Structure-Based Virtual Screening for Drug Discovery: A Problem-Centric Review. The AAPS journal 2012, 14, 133–141. (2) Lavecchia, A.; Giovanni, C. D. Virtual Screening Strategies in Drug Discovery: A Critical Review. Current Medicinal Chemistry 2013, 20, 2839–2860. (3) Gowthaman, R.; Lyskov, S.; Karanicolas, J. DARC 2.0: Improved Docking and Virtual Screening at Protein Interaction Sites. PLoS One 2015, 10, e0131612. (4) Chodera, J. D.; Mobley, D. L.; Shirts, M. R.; Dixon, R. W.; Branson, K.; Pande, V. S. Alchemical Free Energy Methods for Drug Discovery: Progress and Challenges. Current Opinion in Structural Biology 2011, 21, 150–160. (5) Gumbart, J. C.; Roux, B.; Chipot, C. Standard Binding Free Energies from Computer Simulations: What Is the Best Strategy? Journal of Chemical Theory and Computation 2014, 9, 794–802. (6) Kaus, J. W.; Harder, E.; Lin, T.; Abel, R.; McCammon, J. A.; Wang, L. How To Deal with Multiple Binding Poses in Alchemical Relative Protein-Ligand Binding Free 22


Page 22 of 32

Page 23 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Energy Calculations. Journal of Chemical Theory and Computation 2015, 11, 2670– 2679. (7) Wang, L.; Wu, Y.; Deng, Y.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.; Dahlgren, M. K.; Greenwood, J. et al. Accurate and Reliable Prediction of Relative Ligand Binding Potency in Prospective Drug Discovery by way of a Modern Free-Energy Calculation Protocol and Force Field. Journal of the American Chemical Society 2015, 137, 2695–2703. (8) Curcio, R.; Caflisch, A.; Paci, E. Change of the Unbinding Mechanism upon a Mutation: A Molecular Dynamics Study of an Antibody-Hapten Complex. Protein Science 2005, 14, 2499–2514. (9) Aird, A.; Wrachtrup, J.; Schulten, K.; Tietz, C. Possible Pathway for Ubiquinone Shuttling in Rhodospirillum Rubrum Revealed by Molecular Dynamics Simulation. Biophysical Journal 2007, 92, 23–33. (10) Colizzi, F.; Perozzo, R.; Scapozza, L.; Recanatini, M.; Cavalli, A. Single-Molecule Pulling Simulations Can Discern Active from Inactive Enzyme Inhibitors. Journal of the American Chemical Society 2010, 132, 7361–7371. (11) Kingsley, L. J.; Lill, M. A. Including Ligand-Induced Protein Flexibility into Protein Tunnel Prediction. Journal of Computational Chemistry 2014, 35, 1748–1756. (12) Pietrucci, F.; Marinelli, F.; Carloni, P.; Laio, A. Substrate Binding Mechanism of HIV-1 Protease from Explicit-Solvent Atomistic Simulations. Journal of the American Chemical Society 2009, 131, 11811–11818. (13) Limongelli, V.; Bonomi, M.; Marinelli, L.; Gervasio, F. L.; Cavalli, A.; Novellino, E.; Parrinello, M. Molecular Basis of Cyclooxygenase Enzymes (COXs) Selective Inhibition. Proceedings of the National Academy of Sciences 2010, 107, 5411–5416.

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(14) Limongelli, V.; Bonomi, M.; Parrinello, M. Funnel Metadynamics as Accurate Binding Free-Energy Method. Proceedings of the National Academy of Sciences 2013, 110, 6358–6363. (15) Torrie, J. M.; Valleau, J. P. Non-Physical Sampling Distributions in Monte-Carlo FreeEnergy Estimation Umbrella Sampling. Journal of Computational Physics 1977, 23, 187–199. (16) Sun, H.; Tian, S.; Zhou, S.; Li, Y.; Li, D.; Xu, L.; Shen, M.; Pan, P.; Hou, T. Revealing the Favorable Dissociation Pathway of Type II Kinase Inhibitors via Enhanced Sampling Simulations and Two-End-State Calculations. Scientific Reports 2015, 5, 8457. (17) Tiwary, P.; Parrinello, M. From Metadynamics to Dynamics. Physical Review Letters 2013, 111, 230602. (18) Tiwary, P.; Limongelli, V.; Salvalaglio, M.; Parrinello, M. Kinetics of Protein-Ligand Unbinding: Predicting Pathways, Rates, and Rate-Limiting Steps. Proceedings of the National Academy of Sciences 2015, 112, E386–E391. (19) Shan, Y.; Kim, E. T.; Eastwood, M. P.; Dror, R. O.; Seeliger, M. A.; Shaw, D. E. How Does a Drug Molecule Find Its Target Binding Site? Journal of the American Chemical Society 2011, 133, 9181–9183. (20) Buch, I.; Giorgino, T.; De Fabritiis, G. Complete Reconstruction of an EnzymeInhibitor Binding Process by Molecular Dynamics Simulations. Proceedings of the National Academy of Sciences of the United States of America 2011, 108, 10184–10189. (21) Huang, D.; Caflisch, A. The Free Energy Landscape of Small Molecule Unbinding. PLoS Computational Biology 2011, 7, e1002002. (22) Plattner, N.; Noé, F. Protein Conformational Plasticity and Complex Ligand-Binding

24


Page 24 of 32

Page 25 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Kinetics Explored by Atomistic Simulations and Markov Models. Nature Communications 2015, 6, 7653. (23) Doerr, S.; De Fabritiis, G. On-the-Fly Learning and Sampling of Ligand Binding by High-Throughput Molecular Simulations. Journal of Chemical Theory and Computation 2014, 10, 2064–2069. (24) Dickson, A.; Brooks III, C. L. WExplore: Hierarchical Exploration of High-Dimensional Spaces Using the Weighted Ensemble Algorithm. The Journal of Physical Chemistry B 2014, 118, 3532–3542. (25) Huber, G.; Kim, S. Weighted-Ensemble Brownian Dynamics Simulations for Protein Association Reactions. Biophysical Journal 1996, 70, 97–110. (26) Fisher, R. Dispersion on a Sphere. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 1953, 217, 295–305. (27) Sra, S. A Short Note on Parameter Approximation for Von Mises-Fisher Distributions: And a Fast Implementation of Is(x). Computational Statistics 2012, 27, 177–190. (28) Harding, M. W.; Galat, A.; Uehling, D. E.; Schreiber, S. L. A Receptor for the ImmunoSuppressant FK506 Is a Cis-Trans Peptidyl-Prolyl Isomerase. Nature 1989, 341, 758– 760. (29) Siekierka, J. J.; Hung, S. H. Y.; Poe, M.; Lin, C. S.; Sigal, N. H. A Cytosolic Binding Protein for the Immunosuppressant FK506 Has Peptidyl-Prolyl Isomerase Activity but Is Distinct from Cyclophilin. Nature 1989, 341, 755–757. (30) Duyne, G. D. V.; Standaert, R. F.; Karplus, P.; Schreiber, S. L.; Clardy, J. Atomic Structures of the Human Immunophilin FKBP-12 Complexes with FK506 and Rapamycin. Journal of Molecular Biology 1993, 229, 105–124.

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(31) Griffith, J. P.; Kim, J. L.; Kim, E. E.; Sintchak, M. D.; Thomson, J. A.; Fitzgibbon, M. J.; Fleming, M. A.; Caron, P. R.; Hsiao, K.; Navia, M. A. X-Ray Structure of Calcineurin Inhibited by the Immunophilin-Immunosuppressant FKBP12-FK506 Complex. Cell 1995, 82, 507–522. (32) Kissinger, C. R.; Parge, H. E.; Knighton, D. R.; Lewis, C. T.; Pelletier, L. A.; Tempczyk, A.; Kalish, V. J.; Tucker, K. D.; Showalter, R. E.; Moomaw, E. W. et al. Crystal Structures of Human Calcineurin and the Human FKBP12-FK506-Calcineurin Complex. Nature 1995, 378, 641–644. (33) Choi, J.; Chen, J.; Schreiber, S. L.; Clardy, J. Structure of the FKBP12-Rapamycin Complex Interacting with Binding Domain of Human FRAP. Science 1996, 273, 239– 242. (34) Galat, A. Functional Diversity and Pharmacological Profiles of the FKBPs and Their Complexes with Small Natural Ligands. Cellular and Molecular Life Sciences 2013, 70, 3243–3275. (35) Fujitani, H.; Tanida, Y.; Ito, M.; Jayachandran, G.; Snow, C. D.; Shirts, M. R.; Sorin, E. J.; Pande, V. S. Direct Calculation of the Binding Free Energies of FKBP Ligands. Journal of Chemical Physics 2005, 123, 084108. (36) Burkhard, P.; Taylor, P.; Walkinshaw, M. D. X-Ray Structures of Small LigandFKBP Complexes Provide an Estimate for Hydrophobic Interaction Energies. Journal of Molecular Biology 2000, 295, 953–962. (37) Brooks, B.; Brooks III, C.; MacKerell Jr, A.; Nilsson, L.; Petrella, R.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S. et al. CHARMM: The Biomolecular Simulation Program. Journal of Computational Chemistry 2009, 30, 1545–1614. (38) Vanommeslaeghe, K.; MacKerell, A. D. Automation of the CHARMM General Force

26


Page 26 of 32

Page 27 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Field (CGenFF) I: Bond Perception and Atom Typing. Journal of Chemical Information and Modeling 2012, 52, 3144–3154. (39) Vanommeslaeghe, K.; Prabhu Raman, E.; MacKerell, A. D. Automation of the CHARMM General Force Field (CGenFF) II: Assignment of Bonded Parameters and Partial Atomic Charges. Journal of Chemical Information and Modeling 2012, 52, 3155–3168. (40) Eastman, P.; Friedrichs, M. S.; Chodera, J. D.; Radmer, R. J.; Bruns, C. M.; Ku, J. P.; Beauchamp, K. A.; Lane, T. J.; Wang, L.-P.; Shukla, D. et al. OpenMM 4: A Reusable, Extensible, Hardware Independent Library for High Performance Molecular Simulation. Journal of Chemical Theory and Computation 2013, 9, 461–469. (41) Dickson, A.; Mustoe, A. M.; Salmon, L.; Brooks III, C. L. Efficient In Silico Exploration of RNA Interhelical Conformations Using Euler Angles and WExplore. Nucleic Acids Research 2014, 42, 12126–12137. (42) Dickson, A.; Warmflash, A.; Dinner, A. R. Separating Forward and Backward Pathways in Nonequilibrium Umbrella Sampling. The Journal of Chemical Physics 2009, 131, 154104. (43) Dickson, A.; Maienschein-Cline, M.; Tovo-Dwyer, A.; Hammond, J. R.; Dinner, A. R. Flow-Dependent Unfolding and Refolding of an RNA by Nonequilibrium Umbrella Sampling. Journal of Chemical Theory and Computation 2011, 2710–2720. (44) Vanden-Eijnden, E.; Venturoli, M. Exact Rate Calculations by Trajectory Parallelization and Tilting. The Journal of Chemical Physics 2009, 131, 044120. (45) Wood, A. T. Simulation of the Von Mises Fisher Distribution. Communications in Statistics - Simulation and Computation 1994, 23, 157–164.

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 32

(46) Beauchamp, K. A.; Bowman, G. R.; Lane, T. J.; Maibaum, L.; Haque, I. S.; Pande, V. S. MSMBuilder2: Modeling Conformational Dynamics on the Picosecond to Millisecond Scale. Journal of Chemical Theory and Computation 2011, 7, 3412–3419. (47) Senne, M.; Trendelkamp-Schroer, B.; Mey, A. S. J. S.; Sch¨ utte, C.; Noé, F. EMMA: A Software Package for Markov Model Building and Analysis. Journal of Chemical Theory and Computation 2012, 8, 2223–2238. (48) Rojnuckarin, A.; Livesay, D. R.; Subramaniam, S. Bimolecular Reaction Simulation Using Weighted Ensemble Brownian Dynamics and the University of Houston Brownian Dynamics Program. Biophysical Journal 2000, 79, 686–693. (49) Zhang, B. W.; Jasnow, D.; Zuckerman, D. M. Efficient and Verified Simulation of a Path Ensemble for Conformational Change in a United-Residue Model of Calmodulin. Proceedings of the National Academy of Sciences of the United States of America 2007, 104, 18043–18048. (50) Adelman, J. L.; Grabe, M. Simulating Rare Events Using a Weighted Ensemble-Based String Method. The Journal of Chemical Physics 2013, 138, 044105. (51) Saglam, A. S.; Chong, L. T. Highly Efficient Computation of the Basal kon Using Direct Simulation of Protein-Protein Association with Flexible Molecular Models. The Journal of Physical Chemistry B 2015, 120, 117–122. (52) Salentin, S.; Schreiber, S.; Haupt, V. J.; Adasme, M. F.; Schroeder, M. PLIP: Fully Automated Proteinligand Interaction Profiler. Nucleic Acids Research 2015, 43, W443– W447. (53) Hubbard, R. E.; Haider, M. K. Hydrogen Bonds in Proteins: Role and Strength. eLS [Online];

John Wiley & Sons, Ltd:

Chichester, 2010;

10.1002/9780470015902.a0003011.pub2.

28


pp 1-7. DOI:

Page 29 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(54) Gallivan, J. P.; Dougherty, D. A. Cation-Pi Interactions in Structural Biology. Proceedings of the National Academy of Science 1999, 96, 9459–9464. (55) Maragliano, L.; Vanden-Eijnden, E. A Temperature Accelerated Method for Sampling Free Energy and Determining Reaction Pathways in Rare Events Simulations. Chemical Physics Letters 2006, 426, 168–175. (56) Abrams, C. F.; Vanden-Eijnden, E. Large-Scale Conformational Sampling of Proteins Using Temperature-Accelerated Molecular Dynamics. Proceedings of the National Academy of Sciences of the United States of America 2010, 107, 4961–4966. (57) Zheng, L.; Chen, M.; Yang, W. Random Walk in Orthogonal Space to Achieve Efficient Free-Energy Simulation of Complex Systems. Proceedings of the National Academy of Sciences 2008, 105, 20227–20232. (58) Wu, X.; Brooks, B. R. Self-Guided Langevin Dynamics Simulation Method. Chemical Physics Letters 2003, 381, 512–518. (59) Wang, T.; Duan, Y. Ligand Entry and Exit Pathways in the Beta2-Adrenergic Receptor. Journal of Molecular Biology 2009, 392, 1102–1115. (60) Carlsson, P.; Burendahl, S.; Nilsson, L. Unbinding of Retinoic Acid from the Retinoic Acid Receptor by Random Expulsion Molecular Dynamics. Biophysical Journal 2006, 91, 3151–3161. (61) Prinz, J.-H.; Wu, H.; Sarich, M.; Keller, B.; Senne, M.; Held, M.; Chodera, J. D.; Sch¨ utte, C.; Noé, F. Markov Models of Molecular Kinetics: Generation and Validation. The Journal of Chemical Physics 2011, 134, 174105. (62) Rao, F.; Caflisch, A. The Protein Folding Network. Journal of Molecular Biology 2004, 342, 299–306.

29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(63) Dickson, A.; Brooks III, C. L. Native States of Fast-Folding Proteins Are Kinetic Traps. Journal of the American Chemical Society 2013, 135, 4729–4734. (64) Copeland, R. A.; Pompliano, D. L.; Meek, T. D. Drug-Target Residence Time and Its Implications for Lead Optimization. Nature Reviews. Drug Discovery 2006, 5, 730–739. (65) Lu, H.; Tonge, P. J. Drug-Target Residence Time: Critical Information for Lead Optimization. Current Opinion in Chemical Biology 2010, 14, 467–474. (66) Pan, A. C.; Borhani, D. W.; Dror, R. O.; Shaw, D. E. Molecular Determinants of Drug-Receptor Binding Kinetics. Drug Discovery Today 2013, 18, 667–673.

30


Page 30 of 32

Page 31 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry DMSO

BUT

DSS

κ=5.8±0.9

κ=4.0±0.2

κ=5.5±0.9

31


DMSO The Journal of BUT Physical Chemistry PageDSS 32 of 32

κ=5.8±0.9

1 2 3 4 5 6 7

κ=4.0±0.2

κ=5.5±0.9


Ligand Release Pathways Obtained with WExplore: Residence Times

Recommend Documents