Unbinding Kinetics of a p38 MAP Kinase Type II Inhibitor from

Subscriber access provided by HACETTEPE UNIVERSITESI KUTUPHANESI

Article

Unbinding kinetics of a p38 MAP kinase type II inhibitor from metadynamics simulations Rodrigo Casasnovas, Vittorio Limongelli, Pratyush Tiwary, Paolo Carloni, and Michele Parrinello J. Am. Chem. Soc., Just Accepted Manuscript • DOI: 10.1021/jacs.6b12950 • Publication Date (Web): 14 Mar 2017 Downloaded from http://pubs.acs.org on March 14, 2017

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.

Journal of the American Chemical Society 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 22

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

Journal of the American Chemical Society

Unbinding kinetics of a p38 MAP kinase type II inhibitor from metadynamics simulations Rodrigo Casasnovasa, Vittorio Limongelli†b,c, Pratyush Tiwary†d, Paolo Carloni*a, Michele Parrinello*e a

Computational Biomedicine (IAS-5/INM-9), Forschungszentrum Jülich, 52425 Jülich, Germany b Università della Svizzera Italiana (USI), Faculty of Informatics, Institute of Computational Science - Center for Computational Medicine in Cardiology, via G. Buffi 13, CH-6900 Lugano, Switzerland. c Department of Pharmacy, University of Naples “Federico II”, via D. Montesano 49, I-80131 Naples, Italy. d Department of Chemistry, Columbia University, 10027 New York, United States of America e Department of Chemistry and Applied Biosciences, ETH Zurich , and Faculty of Informatics, Institute of Computational Science, Università della Svizzera Italiana, via G. Buffi 13, CH6900 Lugano, Switzerland

1 Environment ACS Paragon Plus


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

Keywords: enhanced sampling, protein-ligand unbinding, kinetics, drug design

Abstract Understanding the structural and energetic requisites of ligand binding towards its molecular target is of paramount relevance in drug design. In the recent years atomistic free energy calculations have proven to be a valid tool to complement experiments in characterizing the thermodynamic and kinetic properties of protein/ligand interaction. Here, we investigate, through a recently developed metadynamics-based protocol, the unbinding mechanism of an inhibitor of the pharmacologically relevant target p38 MAP kinase. We provide a thorough description of the ligand unbinding pathway identifying the most stable binding mode and other thermodynamically relevant poses. From our simulations, we estimated the unbinding rate as koff = 0.020 ± 0.011 s-1. This is in good agreement with the experimental value (koff = 0.14 s-1). Next, we developed a Markov state model that allowed identifying the rate-limiting step of the ligand unbinding process. Our calculations further show that the solvation of the ligand and of the active site play crucial roles in the unbinding process. This study paves the way to investigations on the unbinding dynamics of more complex p38 inhibitors and other pharmacologically relevant inhibitors in general, demonstrating that metadynamics can be a powerful tool in designing new drugs with engineered binding/unbinding kinetics.


Page 2 of 22

Page 3 of 22

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


Introduction Modern drug design strategies attempt at identifying ligands not only with large affinity for their protein targets (measured by the thermodynamic equilibrium binding constant Kb) but also with controlled residence times.1,2 In particular, drug selectivity in vivo may be greatly improved if the residence time for the drug target is large and, more specifically, when it is larger than those of unwanted protein targets.3 In this way, drugs may unbind from undesired proteins for which they have significant affinity and remain on their desired target after some time they are administered.3 However, while highly desirable, manipulating drugs’ chemical structure to tune unbinding kinetics has proven to be very difficult.4,5 While binding rate constants (kon) can be approximated to some extent via a variety of methodsa, the statistically accurate prediction of the pharmaceutically relevant koff along with a structural characterization of the transition state is extremely challenging. A possible path to the calculation of unbinding rate koff through computer simulation is to first estimate the binding rate kon and then relate it to koff using the relation kon / koff = Kb, where Kb is the binding constant. This requires simulating a very large number of binding events to obtain a reliable statistics. In spite of valiant efforts,8,9,12 this has proven challenging. In addition, the error in kon adds to the error in Kb leading to even larger uncertainties. A direct calculation of koff would appear to be a more natural route. Unfortunately, the presence of large free energy barriers leads to very large unbinding times and frustrates any direct calculation of koff. Several methods have been suggested to solve this problem.13-19 For example, Markov State Models are used to extrapolate the long time behavior of the system from a large number of relatively short MD simulations Using this approach, both kon and koff have been calculated.13 Very recently, a different approach based on metadynamics15,16 has been proposed to calculate rates from rare events such as those that are to be found in ligand unbinding.17-19 This approach is similar in spirit to potential flooding20 and hyperdynamics21 and has proven successful in a number of applications including the calculation of koff of benzamidine from the trypsin binding site18 or the study of the water role in the kinetics of cavity-ligand unbinding.22 In this paper, which is meant to be the first of a series, we apply the method to the more pharmaceutically relevant case of mitogen-activated protein kinase p38. This

enzyme is a serine/threonine kinase that controls the biosynthesis of cytokines and is involved in the onset of inflammatory processes.23,24 Designing inhibitors of this kinase represents an attractive therapeutic strategy to counteract inflammation and related diseases. Here, we study the unbinding of 1-(3-(tert-butyl)-1-(p-tolyl)-1Hpyrazol-5-yl)urea (I, Figure 1) that is the basic scaffold of a series of potent inhibitors of p38 whose rates have been experimentally measured.4,25 Our study predicts a koff in agreement with experiments and provides an atomistic picture of the unbinding mechanism, also characterizing the rate-limiting step. From this information other and possibly more potent inhibitors can be designed.



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

Methods Simulation setup. The structure of the human p38 MAP Kinase in complex with I is not available. To model it we started from the structure PDB id 1KV2 of p38 in complex with the structurally similar BIRB-796 molecule (Figure 1) by replacing the 4-(2morpholinoethoxy)napthyl moiety with a hydrogen atom. Since in this crystallographic structure (PDB id 1KV2) the protein lacks two loops involving residues Asn115-Leu122 and Gly170-Ala184, we modeled them by using the Modeller software.26 The p38-I complex was inserted in a box of dimensions 98.4, 81.3 and 75.9 Å, solvated with 14,707 water molecules. In this way each protein atom is at least 24 Å away from its periodic image. 19 Na+ and 12 Cl- ions were added to neutralize the charge of the entire system and to reproduce the ionic strength of the experiments in which the dissociation constants (koff) were measuredb.4,25 A neutral protonation state for ligand I was chosen according to its pKa value27 and to the conditions of the unbinding experiments.4,25 Accordingly to the experimental pH conditions, the aspartate and glutamate residues were considered in their deprotonated state. The lysine and arginine residues were considered in the protonated state and all the other residues with acidbase properties were considered in the neutral state, including all the histidines. The latter were protonated in their ε-nitrogen atom. The final system contained a total of 49,813 atoms. For the protein and ions the AMBER ff99SB-ILDN28-31 force field was used. Water was described by the TIP3P32 model while the ligand force field was constructed using the General Amber Force Field (GAFF).33 The partial atomic charges of I were obtained from HF/6-31G(d) RESP calculations, consistently with the GAFF parametrization.33 All the simulations were performed with the NAMD 2.934 code patched with Plumed2.1.35 Periodic boundary conditions were used. The non-bonding interactions were calculated with a PME scheme with a 10 Å cut-off for the part in real space. All the bonds involving hydrogen atoms were restrained with the SETTLE algorithm.36 The integration time step was set to 2 fs. 1 ns of Langevin dynamics with a friction coefficient of 1 ps-1 was carried out to thermalize the system at 300 K. Subsequently, 5 ns of Langevin dynamics were performed with the Nosé-Hoover barostat37,38 until the average pressure of the system was equilibrated to 1 atm. Because the initial structure of the p38-I complex was modeled from X-ray experiments together with two protein loops, a long equilibration of 500 ns was performed, during which the protein RMSD equilibrated. Metadynamics simulations. Metadynamics16,39,40 is a well-known and much used enhanced sampling method, whose validity has been formally established.41 More recently it has been shown that this method can be employed to obtain reaction rates.17,19 This variant of metadynamics, which falls under the name of infrequent metadynamics will be used here. Like many other enhanced sampling methods, metadynamics is based on the identification of suitable collective variables (CVs), whose fluctuations are gradually enhanced through the course of the simulation, aiding in escape from stable configurations where the system would normally be trapped. The choice of the CVs is very important for a correct simulation. However, as we shall see in our case, the choice of two different sets of CVs led to very similar results, giving added confidence on the reliability of our approach. The set of CVs used were: (i) This set is comprised of two collective variables, CV1 and CV2. CV1 is the distance between the ligand and the binding pocket, defined respectively by the center of mass of the ligand heavy atoms and the center of mass of the alpha carbons of residues His148, Ile146, Ile142, Leu74, Glu71, Leu75, Met78, Ile84, Val83, Phe169, Asp168, Leu167, Ile166, Leu108, Val38, Ala51, Lys53, Thr106,


Page 4 of 22

Page 5 of 22

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


Leu104. Simulations that only used CV1 failed to describe the unbinding kinetics. We traced this to a lack of a proper description of water solvation. Thus, we used a second collective variable CV2 meant to describe the solvation of the ligand. This is defined by the number of water oxygens that are within 8 Å of the ligand center of mass. Notice that CV2 is a simplified definition of the ligand’s solvation. Indeed it does not account for all the waters that are within 8 Å from the ligand’s atoms but from its center of mass. Also, since the CVs have to be defined in terms of smoothly differentiable functions42, the number of water molecules solvating the ligand is calculated with a switching function whose mathematical expression is described in the Supporting Information. (ii) In this second set of simulations we used the path CVs introduced by Branduardi et al.43 and successfully used in ligand/protein binding studies.44-46 These are usually named s and z, where s is the progression along a reference path and z is the distance orthogonal to the path.43 In order to find a reference path, we first performed a steered MD simulation47 on the distance between the ligand and the binding pocket (i.e. CV1). Recently two different pathways have been described for the unbinding of p38 inhibitors.48 In our simulations, steering on CV1 does not presume a particular unbinding pathway for the unbinding. Several SMD simulations with different steering velocities and forces very similar unbinding pathways. From the SMD trajectories we chose 8 frames describing the unbinding of I from the bound pose to the unbound state to build initial path CVs. Then, a preliminary rough metadynamics was performed on these initial path CVs and 8 new frames were chosen from the resulting trajectory to define a better path. Although solvation state was not explicitly introduced in this path, it turned out to play an important role in these runs. We note that the use of path CVs allows visiting states even far from the initial path, thus ensuring the exploration of all the possible ligand unbinding pathways. Notice that these two sets of CVs are significantly uncorrelated with each other. The first set involves an explicit accounting of the solvation state. The second set, namely the path CVs, has no explicit information about solvation. Given that non-trivial and slow fluctuations in solvation have been expected to be one of the molecular determinants of unbinding9,49, these two sets are thus qualitatively very different by construction. Furthermore, in the path CV framework, the inclusion of the z variable facilitates a thorough exploration of the configuration space orthogonal to the initial chosen path, a feature missing in the first set of CVs. Further information on these two sets of CVs can be found in the Supporting Information. Calculation of koff from metadynamics. Taking cue from the work of Grubmüller20 and Voter21, Tiwary and Parrinello noticed that metadynamics, that was originally meant only to calculate static properties, could be used to calculate transition rates in a rare event scenario.19 The crucial assumption for this scheme to work is that zero bias is added in the transition state region, and that the biased CVs can distinguish between all metastable states of relevance. In the work of Tiwary and Parrinello, this goal was obtained by an infrequent deposition of the Gaussians.19 It was then later shown that the assumptions on which the method is based can be checked by working a statistical analysis aimed at establishing if the distribution of escaping times is Poissonian,17 a property that should be satisfied if the rare event hypothesis holds.17 From a metadynamics run one can relate the time measured in metadynamics tMTD to the real one t via the equation: =


[1]


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 6 of 22

where α is an acceleration factor that can be computed from the added metadynamics bias V(s,tMTD): = 〈 ( , ) 〉

[2]

where β is the inverse of kbT and s are the values of the CVs at the time tMTD. In practice, the following equation can be used to calculate the rescaled time t at any given point of the simulation. Hence, the rescaled time t is calculated by summing the time steps taken so far rescaled at each step40: = ∑ d ( ( ), )

[3]

where nMTD is the total number of steps in the metadynamics simulation and tiMTD = idt is the metadynamics time at the i-th step40. A statistical analysis based on the two-sample Kolmogorov-Smirnov (KS) test can quantitatively assess how precisely the above assumptions have been met. Thus, if significant bias got deposited in the TS region even with infrequent biasing, or if there are hidden slow motions at play that the biasing CV does not capture, it would lead to failing the test for time-homogeneous Poisson behavior. This would result in a low pvalue on the KS test, where a value higher than 0.05 is traditionally considered satisfactory for Poisson statistics17. The problem of ligand unbinding is especially suited for this approach. Indeed, numerous investigations over the years have provided evidence in favor of an energy landscape with few high and sharp barriers, resulting in an overall Arrhenius behavior of the rates and relatively short travel times through the barriers.50-53 In particular, the experimental koff for ligand I, BIRB-796 and all the other derivatives were derived by considering a simple competition mechanism, which inherently assumes single exponential kinetics, to fit the experimental data.4,25,54 In practice, given a set of CVs, N independent metadynamics simulations with infrequent bias deposition were started from the ligand in the bound pose and stopped when the ligand reached the unbound state. For each of these metadynamics simulations, the rescaled unbinding time is calculated via equation [3]. Then, the characteristic unbinding time τ is calculated as the average of the various unbinding time observations, and the unbinding rate koff as its reciprocal. An empirical cumulative distribution function (ECDF) is built with the N unbinding times calculated from the simulations and a theoretical cumulative distribution function (TCDF) is built from a large number of times randomly generated according to a cumulative distribution function of a homogeneous Poisson process with characteristic time τ17:

TCDF = 1 − exp "− $ #

[4]

Finally, a two sample KS test is performed to test the null hypothesis that the ECDF and the TCDF show no significant differences and correspond to the same underlying distribution. The null hypothesis is rejected if the p-value resulting from the KS test is lower than the significance level α, which is usually taken as 0.05. In the present case, N=17 when we used CV1 plus CV2, N=10 when we used the pathCVs and N=10 when we used only CV1. Another important parameter that influences the quality of the reconstructed dynamics is the frequency of bias deposition in the metadynamics simulation. We are not aware of an a priori procedure to know how infrequently the bias should be deposited. Instead, we make use of the KS test to a posteriori assess if the frequency of bias deposition was appropriate. In practice, we performed several tests starting with


Page 7 of 22

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


bias deposition of 1 ps-1 and reducing down to 0.1 ps-1, from which we obtained the results reported in this article. Note that in this work, metadynamics with frequent biasing in the conventional sense is not used to obtain the free energy profile and the corresponding thermodynamics of binding. Here instead, metadynamics is used to enhance the fluctuations of the bound ligand and promote its unbinding. Also note that, according to the Tiwary and Parrinello formulation,19 the transition state regions should not be modified with the metadynamics bias in order to allow a correct calculation of the rates. Hence, the infrequent bias deposition is used during our runs to avoid that the transition state regions are corrupted by such bias during the transition between stable states. In addition, we stop the simulations once the ligand reaches the unbound state for the first time to avoid that further recrossings to the bound state corrupt the transition state region. This scheme implies that neither the transition state nor the unbound state regions of the free energy surface are thoroughly sampled and converged. Therefore, this protocol does not provide converged free energies for the calculation of binding thermodynamics. Unbinding mechanism through construction of a Markov model. We have built a kinetic model for the state-to-state transitions to understand in detail the unbinding mechanism. The states of the model were identified using CV1 and CV2 as coordinates from the full unbinding metadynamics calculations. These provide information on where the system spends most of the time18,19. Following ref.18, we considered only those states in which the system spends much longer times than successive bias deposition during the metadynamics run. In the present case, the bias deposition period was 10 ps so, only configurations where the system spends a minimum of MD simulation time of 1 ns or more were considered. Comparable nanosecond lag times have been also used in Markov model studies of protein-ligand association and protein folding. 12,13,55 Apart from the unbound state (US), 3 such configurations were found, labeled as BS (bound state) and intermediates S2, S3. In order to identify any kinetically relevant states that the previous simple classification might have missed due to an imperfect choice of CV, we performed RMSD clustering on the three states identified that the system visits before the unbinding. For this purpose, we aligned the trajectories by using the protein heavy atoms as a reference. Then, we performed a RMSD-based cluster analysis56 with a cutoff of 2 Å on the ligand heavy atoms. Finally, the states were drawn from the most populated clusters where again the system spends a minimum time of 1 ns. All the trajectories of the independent full unbinding metadynamics biasing CV1 and CV2 were used for this analysis. Six states were identified (i.e. BS, S2, S3, aBS, aS2, US). Then, we performed an additional collection of independent metadynamics simulations starting from each of the states with infrequent deposition of bias on CV1 and CV2. Here, each independent metadynamics simulation was stopped when the system escapes from the starting state for the first time. This event can be easily detected monitoring the CVs as function of the simulation time. We used the original 17 independent full unbinding simulations to obtain the escape times from BS. In this case, the times were only counted until the system left BS instead of counting until it reached US. For S2, aBS and aS2, 14, 16 and 14 independent simulations were performed respectively in each case. For S3 we could perform more independent simulations (i.e. 32) because the system escaped in relatively shorter computing time. Note that no such runs were performed for the US since this model was designed keeping in mind absorbing boundary conditions in the US. In this case, the rescaled times of each simulation were calculated again with equation [3] and the lifetime of each state was calculated by averaging the times that the independent simulations took to escape such state. Again, the respective KS tests were run for all each group of independent metadynamics simulations that were started from



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 same state. The individual kinetic constants kij, describing the transition from state i to state j, were calculated by taking into account the lifetime of state i and the number of transitions from this state to state j as in refs.18,57. The state-to-state kij constants were used to build a matrix for the full kinetic process (Table S1). The negative of the largest non-zero eigenvalue of the transpose matrix gives the fundamental relaxation rate of the kinetic model58, which in our case corresponds to an estimate of the koff. The associated eigenvector gives an estimate of the rate-determining step (Figure S1). Results Unbinding of I. As introduced previously we use the variant of metadynamics described in ref.19 and inspired by the work of Grubmüller20 and Voter21. In this approach the metadynamics bias is updated so infrequently so as to leave the transition state free of bias, such that the metadynamics transition rates can be related to the real ones by a simple relation that involves the deposited bias. This multiplicative factor that transforms metadynamics times into rescaled times is called acceleration factor. Furthermore, the statistical quality of the times estimated can be assessed by a parameter p that measures the distance of the distribution of the escape times from the Poissonian distribution that characterizes a rare event.17 In the present case, this metadynamics protocol was applied (see Methods) to study the unbinding of I from kinase p38. Since there is no X-ray structure for the holo form, we had to start from a model of the ligand binding mode discussed in detail in the Methods section. Different sets of collective variables (CVs) were used in independent metadynamics simulations to drive the ligand out of the binding site. We first used the protein-ligand distance as CV1 and the ligand solvation as CV2. These full unbinding calculations led to the estimate of koff = 0.02 ± 0.01 s-1 with a p-value = 0.40 in the KS test (Table 1) for the full unbinding. A similar unbinding rate (koff = 0.04 ± 0.03 s-1, p-value = 0.77, Table 1) was obtained when performing metadynamics using path CVs that describe the contact between the ligand and the protein (see Methods for details). These two independent estimations are in good agreement with the value reported from experiment25 (0.14 s-1), especially considering the uncertainty of force-fields. We carried out an additional full unbinding metadynamics simulations where we used only the protein-ligand distance (CV1). In such case, the koff estimate is equal to (4 ± 4)·10-8 s-1, which is several orders of magnitude different from our previous estimates and the experimental value. The very low p-value in the KS test (p-value 0.02) is indeed an indication that the used CV does not properly take into account all the relevant slow degrees of freedom for ligand unbinding.17 Nevertheless, as mentioned earlier, good estimations of koff were obtained when biasing the path CVs without explicitly introducing the solvation coordinate in the definition of the path (Table 1). This brings out one of the most outstanding features of the metadynamics method – there is no one precise choice of CV that needs to be identified for efficient sampling, as recently shown by Tiwary and Berne59,60. This flexibility in the choice of CVs, and the capability to a posteriori statistically quantify the reliability of the dynamics so reconstructed, should be of much use in sampling the dynamics of a wide variety of complex molecular systems. Table 1. Unbinding constants calculated in metadynamics for different collective variables. Experimental i. Distance, Solvation ii. Path (S, Z) iii. Distance aThe

koff (s-1)a 0.14 0.020 ± 0.011 0.04 ± 0.03 (4 ± 4)·10-8

p-value n.a. 0.40 0.77 0.02

calculated koff values are reported as mean ± standard error of the mean.


Page 8 of 22

Page 9 of 22

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


For simplicity we describe the complex ligand unbinding as a sequence of transitions between states (Figure 2, Figure 3). These are sometimes similar to real poses while in others they are a set of states that share common features. In details, we find that starting from the bound state (BS) the ligand leaves the binding site towards the solvent and reaches a metastable state (S2). Here, assisted by water molecules, the ligand’s urea group rotates assuming a more solvent exposed position. The assistance of water molecules is crucial to weaken the hydrophobic interactions formed between the ligand’s t-butyl group and the hydrophobic residues present in the binding pocket. This event is preliminary for the motion of the ligand towards the fully solvated state (US). In the following paragraphs, we describe state by state the whole unbinding process, providing a detailed analysis of the structures of the energetically stable poses and the transition states connecting the stable and metastable states. Bound state (BS). This is the most stable and long-lived state (Figure 2). The mean lifetime of I in the BS is 51 ± 29 s. As explained in the Methods section, this was calculated from the mean escape rescaled times from BS of 17 independent metadynamics simulations. In the BS, the urea group of I forms H-bond interactions with the backbone amide of Asp168 and the carboxylate sidechain of Glu71 (Figure 3). Additional hydrophobic interactions are engaged by the t-butyl group of I in the pocket formed by the residues Leu74, Leu75, Met78, Val83, Ile84, Ile141, Ile146, His148, Ile166 and Leu167 (Figure 3). The toluene ring of I forms additional hydrophobic interactions with residues in the αC helix of the kinase N-lobe. These interaction points are made with the methylene groups of Arg70 and Glu71 sidechains and to a lesser extent the Arg67 sidechain that shows indeed a high conformational flexibility. The methylene carbons of Asp168 sidechain and the backbone atoms of Phe169 in the activation loop form additional hydrophobic interactions with the toluene group. The unsubstituted nitrogen atom of the pyrazole ring faces the exit of the cavity and does not form any Hbond interaction with the protein but it is in contact with few solvent molecules. From BS to State 2 (TS1). While leaving the bound state, the ligand moves its urea group towards the solvent by pivoting on the t-butyl group. This transition, TS1, leads to a second stable state (S2). In order to better characterize the ensemble of TS1 we performed a committor analysis61 (see Supporting Information) on the metadynamics trajectories that started on BS and lead successfully to S2. In the resulting ensemble, the distance between the centers of mass of the ligand and the binding pocket (CV1) is 7.4 ± 0.3 Å and the number of water molecules solvating the ligand (CV2) is 1.4 ± 0.2. From the structural point of view, in TS1 the H-bonds between I’s urea and Asp168 and Glu71 are no longer a direct interaction, but they are mediated by water molecules (Figure 3). This is due to a pivoting motion of the ligand around the t-butyl group and the hydrophobic pocket. In TS1 these hydrophobic interactions are preserved and no water molecules are found between the t-butyl group and the pocket residues. State 2 (S2). State 2 shows a lifetime of 19 ± 7 microseconds. In this state, the ligand I is tilted towards the exit of the channel allowing the entry of water molecules in the binding cavity. The hydrophobic interactions of the t-butyl group are maintained while the urea group forms H-bonds with water molecules and also is able to engage direct or water-mediated H-bonds with residues at the entrance of the pocket like Lys66, Arg67, Arg70, Glu71, Arg149, Asp168, Gly170, Ala172, His174, Gln325 and Glu328. From S2 to State 3 (TS2). A second transition occurs when the ligand leaves S2 to reach a more solvent exposed state S3. During this motion, the t-butyl group of I leaves the hydrophobic pocket that is rapidly filled by few water molecules (Figure 3). The rolling motion of the t-butyl group on the protein surface allows the progressive hydration of the ligand and the binding pocket. A similar role, played by the water molecules in weakening and disrupting ligand/protein interactions during ligand unbinding has been also found in other systems.22,62-65



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

State 3 (S3). In state 3 the ligand resides at the mouth of the binding cavity, in the region formed by the residues Arg70, Leu74, Asp145, Ile146, Ile147, Arg149, His174, Tyr200, Tyr323, Gln325 and Glu328 (Figure 2). Although these residues account for 90% of the interactions with the ligand, this state consists of an ensemble of interaction poses in which the ligand exchanges rapidly between them. A cluster analysis only on S3 state shows that in the main pose, with a population of 54%, ligand I forms two H-bonds between the urea group and the backbone of Asp145 and hydrophobic interactions with Ile146 and Ile147 (Figure 2). Subsequent poses show similar H-bond and hydrophobic interactions with the aforementioned residues. The binding conformation of the four most populated clusters, accounting for 79% of the total population, is depicted in Figure 2. The residence time for S3, calculated from the multiple state-to-state metadynamics runs, is 0.5 ± 0.4 µs. Unbound state (US). The unbound state was defined, independently of the employed CVs, by the full solvation of the ligand or by the almost full solvation of the ligand where still few H-bond and/or hydrophobic interactions were formed far (i.e. ~30 Å) from the binding pocket with residues at the protein surface. Alternative Bound State (aBS) and alternative State 2 (aS2). In one of the full unbinding metadynamics runs in which CV1 and CV2 were biased, a new bound state was found using the RMSD criterion described in the section ‘Unbinding mechanism through construction of a Markov model.’. This new bound state differs from the starting binding pose based on the X-ray structure of p38 in complex with BIRB-796 (see Methods section) and to distinguish it from BS, this new bound state will be denoted as alternative bound state (aBS). In this run, starting from BS the ligand I reaches state S2 and, due to its relatively high mobility in such state, undergoes a rotation to bind the pocket in an alternative state 2 (aS2). Due to such a rotation, in aS2 the toluene group of I occupies the hydrophobic pocket while the t-butyl group points towards the solvent (Figure 4). Then, from aBS the ligand moves deeper in the binding pocket to reach aBS. The lifetimes, calculated from the state-to-state independent metadynamics runs, are 0.15 ± 0.09 s for aBS and 0.4 ± 0.2 ms for aS2. In aBS, one nitrogen of I’s urea forms an H-bond with the backbone carbonyl of Ile166 while the other nitrogen forms an H-bond with the backbone carbonyl of Asp165. The tolyl aromatic ring forms hydrophobic interactions with the sidechain of Val83 on one side of the ring and with the sidechain of Met78 and Leu75 on the other side. Additional hydrophobic interactions are formed with the sidechain of Thr140 and the backbone of Lys79, His80 and Gly85. Further hydrophobic interactions are established by the pyrazole ring and the t-butyl group with the sidechain of Leu75 and the backbone of Asp168 and Phe169. In aS2, the urea group of I forms a H-bond interaction with the backbone carbonyl of Ile147. The toluene group forms hydrophobic interactions with the sidechains of Leu74, Leu75, Met78, Ile141, Ile146, His148 and the aliphatic part of Asp168 sidechain. On the other hand, the t-butyl group forms hydrophobic interactions with the sidechain of Leu74, the backbone of Glu71 and the aliphatic part of Arg67, Arg70 and Glu71 sidechains. Unbinding mechanism and rate-limiting step. From the Markov State Model built with the additional state-to-state metadynamics simulations (Figure 3, Table S1) we obtained a relaxation rate for the whole unbinding of 0.002 s-1, which is in fair agreement with the koff estimated with the full unbinding metadynamics simulations. The eigenvectors of the transpose kinetic matrix indicate that the global relaxation of this Markov model corresponds to the system escaping the bound state BS (Figure S1). Since from BS the system evolves to S2 (Figure 3), we conclude that the rate-limiting step is the transition between BS and S2. As shown in Figure S1, the second dominant eigenvector corresponds to the system leaving the alternative state aBS. From the stateto-state rate constants it is possible to estimate the thermodynamic equilibrium


Page 10 of 22

Page 11 of 22

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


constant between BS and aBS using the principle of detailed balance (the detailed procedure of this calculation is reported in the Supporting Information). Expressed in terms of free energy, at 300 K, our estimation is that BS is 5.6 ± 1.3 kcal/mol more stable than aBS. This indicates that I will mostly (i.e. >99.9%) bind p38 in the BS pose. From S2 the ligand reaches S3 and can evolve to the unbound state (US), back to S2 or to aS2 (Figure 3). The eigenvector corresponding to the fastest relaxation mode (λ6, Figure S1) shows that from S3 the most probable transition is to the unbound state US, then to S2 and finally to aS2. The respective probabilities of these transitions are 0.75, 0.16 and 0.09. As shown in Figure 3, from both states S2 and aS2, the rate for the transition to S3 is very similar (i.e. 7.4·103 s-1 vs 2.0·103 s-1 respectively). However the transition of aS2 to aBS is slower than that of S2 to BS (i.e. 5.3·102 s-1 vs 4.4·104 s-1 respectively). For this reason, aS2 presents longer lifetime than S2 (Figure 3). A possible explanation for this is that in the aS2-aBS transition the ligand I breaks and forms both H-bonds and hydrophobic interactions to deepen in the pocket (Figure 4). On the other hand, in the S2 to BS transition only the H-bonds interactions suffer significant rearrangement while the hydrophobic interactions between the t-butyl group and the hydrophobic pocket barely change (Figure 2).

Discussion In the present work we have investigated the unbinding kinetics of a potent kinase inhibitor from its binding site. Our metadynamics-based estimation of koff (0.02 ± 0.01 s1) is in good agreement with experiment (0.14 s-1). Our theoretical model is based on the assumption that the ligand unbinding process presents different energy barriers and one rate-limiting step, thus following an Arrhenius behavior. This hypothesis turns out to be a posteriori statistically validated in the KS tests. In particular, the fact that the distribution of times of the full unbinding metadynamics simulations on CV1+CV2 and the path CVs passed the KS test (i.e. p-value>0.05) is an indication that the unbinding of I can be described with a Poisson distribution with a characteristic transition time. This implies that the unbinding process is dominated by a single well-defined barrier and follows an Arrhenius behavior. Therefore, the infrequent metadynamics method is an appropriate approach for the calculation of the unbinding kinetics in this system. On the other hand, the Markov state model derived from state-to-state metadynamics simulations also indicates that the unbinding is dominated by a single barrier since the global relaxation of the entire model is associated with the transition between BS and S2. Kinetic fluorescence experiments showed that the p38 binding and unbinding of I follow single exponential kinetics.4,25 This is not unique for ligand I, the same type of binding and unbinding kinetics was reported for the other p38 inhibitors in the BIRB796 family.4,25 Arrhenius type of protein-ligand binding and unbinding kinetics has been reported for protein-ligand systems other than the one presented in this work. These include computational and experimental investigations: (i) The work of Huang and Caflisch50 showed that the dissociation of six small ligands from the FK506 binding protein follows a single exponential kinetics because of the presence of a single dominating free energy barrier. (ii) The review of Lee et al.51, reports a variety of experimental studies of protein-ligand unbinding processes that obey an Arrhenius behavior. (iii) Fluorescence competitive dissociation assays and atomic force microscopy experiments indicated single exponential kinetics for ligands binding to nine different anti-fluorescein antibody mutants.52 (iv) Freitag et al.53 performed radiometric competition assays to determine the biotin koff from a streptavidin mutant in a range of temperatures and found a good fit to a single exponential model. (v) Our previous metadynamics study of the unbinding of



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

benzamidine from trypsin also showed an unbinding mechanism comprised by a series of sharp barriers which obeys a single exponential kinetics. The found binding pose aBS shows a significant lifetime compared to that of BS (i.e. 0.15 ± 0.09 s vs 51 ± 29 s). Also, the intermediate aS2 shows larger lifetime than that of S2. Our estimation is that aS2 is slightly more stable than S2 (i.e. 0.5±0.7 kcal/mol). Therefore, it can be that the most probable unbinding mechanism implies the alternative states, that is, they are kinetically relevant for the global koff. However, the koff estimated from the Markov model (0.002 s-1), which explicitly takes into account the alternative states, is in fair agreement with that estimated from the full unbinding metadynamics (0.020 ± 0.011 s-1), where most trajectories did not visit the alternative states. This suggests that the most probable unbinding mechanism proceeds only via intermediates S2 and S3 and does not imply aS2 and aBS. There are further indications supporting this mechanism. First, the pose BS is 5.6 ± 1.3 kcal/mol more stable than the pose aBS so, under thermodynamic equilibrium conditions, the population of I bound in the aBS pose will be less than 1/2000. Second, once I leaves BS and reaches S2, the rate to reach the unbound state US is significantly larger than the rate to reach aS2 (i.e. 2.2·104 s-1 vs 7.7·102 s-1 calculated by taking into account the steady state approximation). Third, which is related to the previous point, is that once I reaches S3, the most probable transition is to US, then to S2 and finally to aS2 with probabilities of 0.75, 0.16 and 0.09 respectively. Consequently, we propose a scheme like BS→S2→S3→US for the most probable unbinding mechanism. Obviously, only an infinite sampling of trajectories ensures that all the unbinding pathways are correctly sampled. In a real case situation with limited number of trajectories, most will sample the minimum free energy path and only a minor fraction of trajectories will sample alternative mechanisms of higher energy. However, missing higher energy pathways in the sampling is not critical since their weight in describing the event decreases exponentially with the energy according to a Boltzmann factor. Our simulations provide a mechanistic explanation of the structure-kinetic relationships reported by Regan et al. from stopped-flow fluorescence experiments for BIRB-796 and its derivatives II-IV.4,25 For convenience, the structure-kinetic data of refs.4,25 is reported in Table 2 for ligands II-IV. Ligand I binds p38 with micromolar affinity, much more weakly than BIRB-796 (i.e. Kd~1.2µM and ~0.1nM respectively).25 Whereas the binding rates for I and BIRB-796 are very similar, the unbinding rate of I is ~105-fold faster. However, the experimental unbinding rates reported by Regan et al. indicate that ligand I is the “minimal” structure for all compounds in the BIRB-796 family to show measurable bound residence times.25 While we have not yet simulated the other ligands II-IV in Table 2, we can speculate from the results obtained in the smallest ligand I about the behavior of the bulkier ligands. Experiments show that the addition of a phenyl (II) or naphthyl (III) groups to the urea nitrogen of I decreases the koff by 102 and 103-fold respectively25 (Table 2). We argue that this may be due to the combined effect of extra stabilization energy due to protein-ligand interactions and of the reduced access of water in the binding cavity that is occupied by larger hydrophobic groups. These are likely to reduce the solvent fluctuations responsible for disrupting the direct protein-ligand interactions and, therefore, increase the lifetime of the bound state. In ligand IV the urea is linked to a pchlorophenyl group and the tolyl group of the pyrazole ring is replaced by a methyl group. Experiments show that ligands I and IV have similar binding affinities for p38 and identical unbinding rates4,25 (Table 2). We hypothesize that the reason of the identical unbinding rates of I and IV is that while the p-chlorophenyl shields the ureaGlu71 H-bond, replacement of the tolyl group by a much smaller methyl group has the opposite effect on the solvent exposed side of the cavity. The ligand binding site in kinases is located close to the activation loop. Kinase type II inhibitors require that a highly conserved Asp-Phe-Gly sequence, known as the DFG triad, of the activation loop undergoes a rearrangement from the active DFG-in


Page 12 of 22

Page 13 of 22

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


conformation capable of binding ATP and type I inhibitors to the inactive DFG-out conformation. 4,5,66 This conformational change creates the hydrophobic pocket that I, BIRB-796 and its derivatives bind in p38 MAP kinase 4,25,66 and other type II inhibitors bind in other kinases.5 NMR and fluorescence kinetic experiments showed that the DFGin to DFG-out transition in p38 MAP kinase occurs in the micro to milliseconds timescale and rate-limits the binding of BIRB-796 and its derivatives.4,66 In our simulations no conformational rearrangements of the DFG motif and the rest of the activation loop are observed for the unbinding of I in any of our metadynamics simulations. Since 1) our estimation of koff of I is in fair agreement with the reported from experiments and 2) if the biased CVs in the metadynamics simulations miss the important slow degrees of freedom, the distribution of unbinding times would fail the KS test, 17 the present simulations suggest that the DFG-out to DFG-in transition is not coupled with the unbinding of I. We would also like to highlight the extensive sampling performed in this work with different classes of biasing collective variables, which is not very commonly done for a system of such complexity. The results reported here represent a total of ~6.8 µs of production simulations among all the full unbinding and the state-to-state metadynamics runs and a total CPU time of 2.5 million core-h were used for this work. The qualitatively and quantitatively consistent nature of the conclusions further increase confidence in the use of metadynamics for studying the unbinding dynamics of such systems. We have shown that metadynamics simulations on the unbinding of I provide valuable insight into the unbinding kinetics and mechanism of the BIRB-796 family of compounds. For this reason, we believe that this technique has great potential to obtain structure-kinetic relationships of other families of inhibitors and help in future drug development based on the optimization of residence times.



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 14 of 22

Table 2. Kinetic and thermodynamic binding/unbinding parameters of I, BIRB-796 and some derivatives. Adapted from refs.4,25. R1

NH

O

NH N N

R2

kon (M-1s-1)

koff (s-1)

Kd (nM)

1.18105

1.410-1

1190

II

7.33104

1.610-3

22

III

1.13105

1.210-4

1.1

1.2105

1.410-1

1160

8.49104

8.310-6

0.098

Ligand I

R1

R2

H

Cl

CH3

IV

BIRB-796

O

N O


Page 15 of 22

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


Figure 1. Left: Chemical structure of ligand I (1-(3-(tert-butyl)-1-(p-tolyl)-1H-pyrazol-5yl)urea) and experimental kon and koff values.4,25 Middle: BIRB-796 and derivatives bind into a cavity formed in the interface of p38 MAP kinase’s N- and C-lobes. Right: ProteinBIRB-796 interactions in the bound state.



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

Figure 2. States BS, S2, S3, US and transition states TS1, TS2 for the unbinding of I. Bulk water is not shown for clarity purposes. Water oxygen atoms are depicted as red spheres. Hydrogen bond interactions are depicted as dashed lines. The centers of mass of the binding pocket and the ligand are depicted as green spheres and are also connected by a dashed line to indicate the distance between the two (i.e. CV1). In the depicted snapshots, the distance between the centers of mass for BS, TS1, S2, TS2 and


Page 16 of 22

Page 17 of 22

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


US are respectively 3.9 Å, 6.9 Å, 8.4 Å, 14.7 Å and 35.4 Å. In S3, from left to right and top to bottom, the distances between centers of mass are in the range 15-20 Å. Note that the center of mass employed for CV1 is related to the entire binding cavity and not only to the residues forming the hydrophobic pocket, which here are depicted with an isosurface.

Figure 3. Representation of the Markov State Model, including the lifetimes of each state and the rates between states computed from the state-to-state metadynamics simulations. All the rates are in s-1.

Figure 4. Alternative bound state (aBS) and alternative state 2 (aS2) for the unbinding of I in which the tolyl group of I occupies the hydrophobic cavity deep in the binding pocket. The states BS and S2 are also shown for the sake of comparison. Water oxygen atoms are depicted as red spheres. Hydrogen bond interactions are depicted as dashed lines. The centers of mass of the binding pocket and the ligand are depicted as green spheres. The centers of mass are also united by a dashed line to indicate the distance between the two (i.e. CV1). In the depicted snapshots, the distance between the centers of mass for aBS and aS2 are respectively 3.1 Å and 8.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

Supporting Information Details on metadynamics simulations, Details on the Markov Model, Tables S1-S4, Figures S1–S3. Corresponding Author information *Corresponding Authors: Michele Parrinello ETH Zurich, Department of Chemistry and Applied Biosciences Via Giuseppe Buffi 13 Lugano, CH 6900 Telephone: +41-58666 48 01 Email: [email protected] Paolo Carloni Forschungszentrum Julich, Computational Biomedicine section, Institute of Advanced Simulation (IAS-5) and at the Institute of Neuroscience and Medicine (INM-9) Wilhelm-Johnen-Straße Julich, Nordrhein-Westfalen, DE 52428 Telephone: +49 2461 61 8941 Email: [email protected] Author contributions †These

authors contributed equally to this work.

Acknowledgements The authors gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JUROPA (Project id HGR26) and JURECA67 (Project id JIAS51) at the Jülich Supercomputing Centre (JSC), Jülich, Germany. VL acknowledges the support from the Swiss National Supercomputing Center (CSCS under project ID s557), the Swiss National Science Foundation (Project N. 200021_163281) and the COST action CA15135 (Multi-target paradigm for innovative ligand identification in the drug discovery process MuTaLig). This project has received funding from the European Union’s Horizon 2020 Framework Programme for Research and Innovation under Grant Agreement No. 720270 (Human Brain Project SGA1). R.C.P. is a recipient of the Human Brain Project SGA1 grant.


Page 18 of 22

Page 19 of 22

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


References

(1) Copeland, R. A. Nat. Rev. Drug Discov. 2016, 15, 87. (2) Copeland, R. A.; Pompliano, D. L.; Meek, T. D. Nat. Rev. Drug Discov. 2006, 5, 730. (3) Walkup, G. K.; You, Z.; Ross, P. L.; Allen, E. K. H.; Daryaee, F.; Hale, M. R.; O'Donnell, J.; Ehmann, D. E.; Schuck, V. J. A.; Buurman, E. T.; Choy, A. L.; Hajec, L.; Murphy-Benenato, K.; Marone, V.; Patey, S. A.; Grosser, L. A.; Johnstone, M.; Walker, S. G.; Tonge, P. J.; Fisher, S. L. Nat. Chem. Biol. 2015, 11, 416. (4) Pargellis, C.; Tong, L.; Churchill, L.; Cirillo, P. F.; Gilmore, T.; Graham, A. G.; Grob, P. M.; Hickey, E. R.; Moss, N.; Pav, S.; Regan, J. Nat. Struct. Mol. Biol. 2002, 9, 268. (5) Schneider, E. V.; Böttcher, J.; Huber, R.; Maskos, K.; Neumann, L. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 8081. (6) Mereghetti, P.; Gabdoulline, R. R.; Wade, R. C. Biophys. J. 2010, 99, 3782. (7) Sandikci, A.; Gloge, F.; Martinez, M.; Mayer, M. P.; Wade, R.; Bukau, B.; Kramer, G. Nat. Struct. Mol. Biol. 2013, 20, 843. (8) Dror, R. O.; Pan, A. C.; Arlow, D. H.; Borhani, D. W.; Maragakis, P.; Shan, Y.; Xu, H.; Shaw, D. E. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 13118. (9) Shan, Y.; Kim, E. T.; Eastwood, M. P.; Dror, R. O.; Seeliger, M. A.; Shaw, D. E. J. Am. Chem. Soc. 2011, 133, 9181. (10) Lane, T. J.; Bowman, G. R.; Beauchamp, K.; Voelz, V. A.; Pande, V. S. J. Am. Chem. Soc. 2011, 133, 18413. (11) Shukla, D.; Meng, Y.; Roux, B.; Pande, V. S. Nat. Commun. 2014, 5. (12) Buch, I.; Giorgino, T.; De Fabritiis, G. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 10184. (13) Plattner, N.; Noe, F. Nat. Commun. 2015, 6. (14) Teo, I.; Mayne, C. G.; Schulten, K.; Lelièvre, T. J. Chem. Theory Comput. 2016, 12, 2983. (15) Barducci, A.; Bussi, G.; Parrinello, M. Phys. Rev. Lett. 2008, 100, 020603. (16) Laio, A.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562. (17) Salvalaglio, M.; Tiwary, P.; Parrinello, M. J. Chem. Theory Comput. 2014, 10, 1420. (18) Tiwary, P.; Limongelli, V.; Salvalaglio, M.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, E386. (19) Tiwary, P.; Parrinello, M. P. Phys. Rev. Lett. 2013, 111, 230602. (20) Grubmüller, H. Phys. Rev. E 1995, 52, 2893. (21) Voter, A. F. Phys. Rev. Lett. 1997, 78, 3908. (22) Tiwary, P.; Mondal, J.; Morrone, J. A.; Berne, B. J. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 12015. (23) Han, J.; Lee, J.; Bibbs, L.; Ulevitch, R. Science 1994, 265, 808. (24) Lee, J. C.; Laydon, J. T.; McDonnell, P. C.; Gallagher, T. F.; Kumar, S.; Green, D.; McNulty, D.; Blumenthal, M. J.; Keys, J. R.; Land vatter, S. W.; Strickler, J. E.; McLaughlin, M. M.; Siemens, I. R.; Fisher, S. M.; Livi, G. P.; White, J. R.; Adams, J. L.; Young, P. R. Nature 1994, 372, 739.



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

(25) Regan, J.; Pargellis, C. A.; Cirillo, P. F.; Gilmore, T.; Hickey, E. R.; Peet, G. W.; Proto, A.; Swinamer, A.; Moss, N. Bioorg. Med. Chem. Lett. 2003, 13, 3101. (26) Šali, A.; Blundell, T. L. J. Mol. Biol. 1993, 234, 779. (27) Barnes, M. J.; Conroy, R.; Miller, D. J.; Mills, J. S.; Montana, J. G.; Pooni, P. K.; Showell, G. A.; Walsh, L. M.; Warneck, J. B. H. Bioorg. Med. Chem. Lett. 2007, 17, 354. (28) Craft, J. W.; Legge, G. B. J. Biomol. NMR 2005, 33, 15. (29) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Proteins: Struct., Funct., Bioinf. 2010, 78, 1950. (30) Pavelites, J. J.; Gao, J.; Bash, P. A.; Mackerell, A. D. J. Comput. Chem. 1997, 18, 221. (31) Walker, R. C.; de Souza, M. M.; Mercer, I. P.; Gould, I. R.; Klug, D. R. J. Phys. Chem. B 2002, 106, 11658. (32) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (33) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157. (34) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. J. Comput. Chem. 2005, 26, 1781. (35) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. Comput. Phys. Commun. 2009, 180, 1961. (36) Miyamoto, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 952. (37) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. J. Chem. Phys. 1995, 103, 4613. (38) Martyna, G. J.; Tobias, D. J.; Klein, M. L. J. Chem. Phys. 1994, 101, 4177. (39) Barducci, A.; Bonomi, M.; Parrinello, M. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 826. (40) Valsson, O.; Tiwary, P.; Parrinello, M. Annu. Rev. Phys. Chem. 2016, 67, 159. (41) Dama, J. F.; Parrinello, M.; Voth, G. A. Phys. Rev. Lett. 2014, 112, 240602. (42) Granata, D.; Camilloni, C.; Vendruscolo, M.; Laio, A. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 6817. (43) Branduardi, D.; Gervasio, F. L.; Parrinello, M. J. Chem. Phys. 2007, 126, 054103. (44) Limongelli, V.; Bonomi, M.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 6358. (45) Limongelli, V.; Bonomi, M.; Marinelli, L.; Gervasio, F. L.; Cavalli, A.; Novellino, E.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 5411. (46) Grazioso, G.; Limongelli, V.; Branduardi, D.; Novellino, E.; De Micheli, C.; Cavalli, A.; Parrinello, M. J. Am. Chem. Soc. 2012, 134, 453. (47) Grubmüller, H.; Heymann, B.; Tavan, P. Science 1996, 271, 997. (48) Sun, H.; Tian, S.; Zhou, S.; Li, Y.; Li, D.; Xu, L.; Shen, M.; Pan, P.; Hou, T. Sci. Rep. 2015, 5, 8457.


Page 20 of 22

Page 21 of 22

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


(49) Pan, A. C.; Borhani, D. W.; Dror, R. O.; Shaw, D. E. Drug Discov. Today 2013, 18, 667. (50) Huang, D.; Caflisch, A. PLoS Comput. Biol. 2011, 7, e1002002. (51) Lee, C.-K.; Wang, Y.-M.; Huang, L.-S.; Lin, S. Micron 2007, 38, 446. (52) Schwesinger, F.; Ros, R.; Strunz, T.; Anselmetti, D.; Güntherodt, H.J.; Honegger, A.; Jermutus, L.; Tiefenauer, L.; Plückthun, A. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 9972. (53) Freitag, S.; Chu, V.; Penzotti, J. E.; Klumb, L. A.; To, R.; Hyre, D.; Le Trong, I.; Lybrand, T. P.; Stenkamp, R. E.; Stayton, P. S. Proc. Natl. Acad. Sci. U. S. A. 1999, 96, 8384. (54) Morelock, M. M.; Pargellis, C. A.; Graham, E. T.; Lamarre, D.; Jung, G. J. Med. Chem. 1995, 38, 1751. (55) Pande, V. S.; Beauchamp, K.; Bowman, G. R. Methods 2010, 52, 99. (56) Daura, X.; Gademann, K.; Jaun, B.; Seebach, D.; van Gunsteren, W. F.; Mark, A. E. Angew. Chem., Int. Ed. 1999, 38, 236. (57) Du, W.-N.; Marino, K. A.; Bolhuis, P. G. J. Chem. Phys. 2011, 135, 145102. (58) Widom, B. J. Chem. Phys. 1971, 55, 44. (59) Tiwary, P.; Berne, B. J. J. J. Chem. Phys. 2016, 145, 054113. (60) Tiwary, P.; Berne, B. J. P. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, 2839. (61) Bolhuis, P. G.; Dellago, C.; Chandler, D. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 5877. (62) Limongelli, V.; Marinelli, L.; Cosconati, S.; La Motta, C.; Sartini, S.; Mugnaini, L.; Da Settimo, F.; Novellino, E.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 1467. (63) Mondal, J.; Friesner, R. A.; Berne, B. J. J. Chem. Theory Comput. 2014, 10, 5696. (64) Mondal, J.; Morrone, J. A.; Berne, B. J. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 13277. (65) Schmidtke, P.; Luque, F. J.; Murray, J. B.; Barril, X. J. Am. Chem. Soc. 2011, 133, 18903. (66) Vogtherr, M.; Saxena, K.; Hoelder, S.; Grimme, S.; Betz, M.; Schieborr, U.; Pescatore, B.; Robin, M.; Delarbre, L.; Langer, T.; Wendt, K. U.; Schwalbe, H. Angew. Chem., Int. Ed. 2006, 45, 993. (67) Krause, D.; Thörnig, P., JURECA: General-Purpose Supercomputer JLRSF 2016, 2, A62 at Jülich Supercomputing Centre, http://dx.doi.org/10.17815/jlsrf-2-121. Footnotes aThese

methods include: Brownian Dynamics (BD) simulations,6,7 Molecular Dynamics (MD) simulations running on optimized hardware,8,9 MD simulations via Markov State Models.10,11 bThe experimental kinetic measurements of refs.4,25 were performed in a buffer containing 20nM Bis-Tris Propane, pH 7.0, 2 mM EDTA, 0.01% (w/v) NaN3 and 0.15% (w/v) n-octylglucoside. The equivalent ionic strength (ca. 34 mM) of NaCl was used in the simulations.



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

TOC Image


Page 22 of 22

Unbinding Kinetics of a p38 MAP Kinase Type II Inhibitor from

Recommend Documents