Enhanced sampling simulations of ligand unbinding kinetics

Subscriber access provided by Macquarie University

Computational Biochemistry

Enhanced sampling simulations of ligand unbinding kinetics controlled by protein conformational changes Yang Zhou, Rongfeng Zou, Guanglin Kuang, Bengt Langstrom, Christer Halldin, Hans Agren, and Yaoquan Tu J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.9b00523 • Publication Date (Web): 27 Aug 2019 Downloaded from pubs.acs.org on August 30, 2019

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 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 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.

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 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

Journal of Chemical Information and Modeling

Enhanced Sampling Simulations of Ligand Unbinding Kinetics Controlled by Protein Conformational Changes

Yang Zhou,a Rongfeng Zou,a Guanglin Kuang,a Bengt Långström,b Christer Halldin, c Hans Ågren, d,a* and Yaoquan Tu a * a

Department of Theoretical Chemistry and Biology, KTH Royal Institute of Technology, AlbaNova University Center, Stockholm, 10691, Sweden b Department of Chemistry, Uppsala University, Uppsala, 75123, Sweden c Department of Clinical Neuroscience, Center for Psychiatric Research, Karolinska Institutet and Stockholm County Council, Stockholm, 17176, Sweden d College of Chemistry and Chemical Engineering, Henan University, Kaifeng, Henan, 475004, P. R. China * Corresponding author: Yaoquan Tu, Department of Theoretical Chemistry and Biology, KTH Royal Institute of Technology, AlbaNova University Center, Stockholm, 10691, Sweden Tel.: +46 8 790 96 45 Email: Y. Tu ([email protected]) Hans Ågren, Department of Theoretical Chemistry and Biology, KTH Royal Institute of Technology, AlbaNova University Center, Stockholm, 10691, Sweden Tel.: +46 8 790 99 83 Email: H. Ågren ([email protected])

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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 Understanding unbinding kinetics of protein-ligand systems is of great importance for the design of ligands with desired specificity and safety. In recent years, enhanced sampling techniques have emerged as effective tools for studying unbinding kinetics of protein-ligand systems at the atomistic level. However, in many protein-ligand systems, the ligand unbinding processes are strongly coupled to protein conformational changes and how to disclose the hidden degrees of freedom closely related to the protein conformational changes so that sampling is enhanced over these degrees of freedom remains a great challenge. Here, we show how potential-scaled molecular dynamics (sMD) and infrequent metadynamics (InMetaD) simulation techniques can be combined to successfully reveal the unbinding mechanism of ASEM (3-(1,4-diazabicyclo[3.2.2]nonan-4-yl)6-[18F]fluo-rodibenzo[b,d]thiophene 5,5-dioxide) from a chimera structure of the α7-nicotinic acetylcholine receptor. By using sMD simulations, we disclosed that the “close” to “open” conformational change of loop C plays a key role in the ASEM unbinding process. By carrying out InMetaD simulations with this conformational change taken into account as an additional collective variable, we further captured the key states in the unbinding process and clarified the unbinding mechanism of ASEM from the protein. Our work indicates that combining sMD and InMetaD simulation techniques can be an effective approach for revealing the unbinding mechanism of a protein-ligand system where protein conformational changes control the unbinding process.


Page 2 of 19

Page 3 of 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


Introduction Rational design of ligands as drugs or detection tracers has traditionally progressed in terms of binding affinities. However, increasingly abundant studies indicate that kinetic properties, such as residence time, binding and unbinding rates, can be as important as, or even more important than, the binding affinity of a ligand for use as a drug or tracer targeting living organisms 1–5. Although the kinetics parameters for ligands can often be measured experimentally 5-6, the way they leave the protein receptor and the molecular determinants that dictate the dissociation process, such as key interactions and conformational changes, still remain largely unknown and can hardly be revealed solely from experiment. Such knowledge is though essential for rational structural optimization for the purpose of improving kinetic properties and further development of new drugs or tracers with unprecedented performance and functionality. Computer simulation has become a powerful tool for studying protein-ligand systems with atomistic resolution. Enabled by the development of methods such as Brownian dynamics, Markov state modeling, and enhanced sampling molecular dynamics (MD) simulation techniques, computer simulation has become progressively more successful in predicting kinetics of protein-ligand systems and in providing atomic-level insight into the binding and unbinding mechanisms of the systems 7–10. Among the methods for studying protein-ligand unbinding kinetics, an enhanced sampling method called infrequent metadynamics (InMetaD) is very efficient in calculating the “koff” rate for a ligand unbinding from its protein receptor. InMetaD has been successfully used for calculating the “koff” rates for systems such as trypsin-benzamidine, biotin-streptavidin, T4 lysozyme-benzene inhibitors from the p38 MAP kinase, and an anticancer drug from the c-Src kinase 11–14. It is possible to use InMetaD to study unbinding processes with residence times up to seconds 12. However, the results from these simulations depend very much on the set of collective variables (CVs) selected. In many protein-ligand systems, the ligand unbinding kinetics are controlled by the protein conformational changes. How to select the CVs closely related to the protein conformational changes so that reliable kinetic properties can be obtained from the corresponding InMetaD simulations remains a great challenge. For the purpose of using InMetaD simulations to study ligand unbinding kinetics, we propose a strategy in which another enhanced sampling method, called potential-scaled MD (sMD), is introduced to disclose the CVs most relevant to the unbinding process. sMD is CV-free. In an sMD simulation, the potential of the system is uniformly scaled down by a scaling factor β 15. This flattens the potential energy surface so that the unbinding process is greatly accelerated. sMD has previously been successfully applied in ranking the koff values of a series of compounds for drug discovery 15,16. Since the koff values are closely related to the unbinding mechanism, the agreement between the sMD results and experiment in ranking the ligands reflects that the unbinding mechanism stays essentially unchanged in the sMD simulations. Indeed, a recent study shows that sMD can be used to reveal a drug unbinding pathway 17. Although the occurrence of unbinding events increases in the sMD simulations, the essential intermediate states during the unbinding process can still be sampled with a relatively high probability. Therefore, we believe that sMD can provide a prior knowledge of the unbinding process and disclose the hidden slow degrees of freedom closely related



to the protein structure changes critical to the process. Here we show how sMD and InMetaD simulation techniques can be combined to reveal the unbinding mechanism of [18F]ASEM (3-(1,4-diazabicyclo[3.2.2]nonan-4-yl)-6-[18F]fluorodibenzo[b,d]thiophene 5,5-dioxide) from a chimera structure of the α7 nicotinic acetylcholine receptor (α7-nAChR)6. The evidence that α7-nAChRs are heavily implicated in a variety of neurological disorders such as schizophrenia, Alzheimer's disease, and traumatic brain injury, has served as a motivation to use positron emission tomography (PET) to image and quantify α7nAChRs in the human brain for early detection of the disorders 18. [18F]ASEM has been shown to have an extraordinarily high binding affinity to the α7-nAChR and is therefore believed to be the most promising PET tracer for in vivo quantification of α7-nAChR in the human central nervous system. However, it has been gradually realized that its slow wash-out from the brain makes the scanning time for accurately quantifying [18F]ASEM too long in humans 6,18. This implies that the “koff” rate for ASEM unbinding from the α7-nAChR could also a critical factor in determining its applicability as a PET tracer. Therefore, understanding the unbinding mechanism of ASEM from the α7-nAChR, and most importantly, identifying what are the critical molecular interactions between ASEM and the α7-nAChR during the unbinding process will help us design PET tracers with desired kinetic properties. In our study, we combined sMD and InMetaD simulation techniques to reveal the unbinding mechanism of ASEM from an α7-nAChR chimera structure named as acetylcholine binding protein (AChBP). The ligand-binding core and flanking regions of the AChBP are lined entirely by the α7-nAChR residues, making this structure a good model for the ligandbinding domain of the α7-nAChR. Because ASEM has an extraordinarily high binding affinity to the α7-nAChR, and in particular, the binding pocket of the AChBP is buried deeply in the protein and covered by a bunch of residences, a study of the unbinding kinetics of ASEM from the AChBP can also serve as a good example of how sMD and InMetaD simulation techniques can be combined to reveal the unbinding mechanism of a ligand from its receptor.

Results Initial infrequent metadynamics simulations In our initial study of the unbinding process of ASEM from the AChBP, two CVs, which have been used in the studies of some drug-target unbinding processes successfully 13,14, i.e. CVdis - the distance between the ligand and the protein pocket, and CVsolv - the number of water molecules in the nearest neighbor solvation shell, were selected (see Method section for the detailed definition of the CVs). Twenty independent InMetaD simulations were carried out, and the residence time from each simulation was rescaled by the acceleration factor calculated according to eq. (4). The Kolmogorov-Smirnov (KS) test was then performed to assess the residence times thus obtained. The p-value from the KS test is 0.003, which is much lower than the confidence threshold 0.05, indicating that the residence times are of a non-Poissonian distribution form. Indeed, as shown in Figure 1a, the obtained empirical cumulative distribution function (ECDF) deviates significantly from the theoretical cumulative distribution function (TCDF) for a Poisson process. According to previous studies13,19, there could be two reasons for that the residence times show a non-Poissonian


Page 4 of 19

Page 5 of 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


distribution: (1) the bias potential is deposited in the transition state region, and (2) some slow degrees of freedom closely related to the unbinding process are not considered. Since the bias potential in our simulations was deposited with a large time lag (10 ps), which is in line with those in other InMetaD simulations 13,14,20, we believe that it is unlikely that the bias potential was deposited in the transition state region. Therefore, we attribute the non-Poissonian distribution to that some CVs closely related to the ASEM unbinding process were not selected, though the distance and solvation CVs have been suggested to play prominent roles in the ligand unbinding processes in several studies 13,14,21. Indeed, Tiwary et al. found that without considering one slow degree of freedom for ligand unbinding, the p-value in the KS test was very low (0.02) and the koff was far from the experimental value 13. After an additional CV was selected to take into account of another hidden slow degree of freedom for the ligand unbinding, the p-value increased, and the estimated koff became in agreement with the experimental value. In our case, since the binding pocket is buried deeply in the protein and covered by a bunch of residues, the hidden slow degrees of freedom could be related to the protein structure changes critical to the ASEM unbinding process, which are, however, very difficult to identify without prior knowledge of the unbinding process.

Figure 1. Cumulative distribution of the residence times from the first round of InMetaD simulations and the representative structures from the InMetaD simulations, sMD simulations and PDB. (a) cumulative distribution of the residence times from 20 independent InMetaD simulations using the distance and solvation CVs. (b) representative structures from the InMetaD simulations: time evolutions of ASEM and loop C are colored from red (t=0) to blue (unbound). (c) representative



configurations in sMD simulations: time evolutions of ASEM and loop C are colored from red (t=0) to blue (unbound). (d) Alignment of representative AChBP structures (PDB: 2BR7, 2C9T, 2XYT, 2Y54, 5XGL). Loop C is colored from red to blue (“close” to “open”). Dashed line is the Cα- Cα distance between TrpB and Y188.

Opening and closing of loop C in the unbinding process disclosed by potential-scaled MD simulations In order to find additional CVs essential for InMetaD simulations of the unbinding process of ASEM, we resort to CV-free sMD simulation to disclose the slow degrees of freedom closely related to the unbinding process. Here, we carried out twenty independent sMD simulations. Comparing the trajectories from sMD and InMetaD simulations, we found that there is a significant difference in the behavior of loop C (Arg182 to Tyr191) that covers the protein binding pocket. As shown in Figure 1c, in the sMD simulations loop C underwent a “close” to “open” conformational change during the ASEM unbinding process, but this conformational change was not observed in the InMetaD simulations (Figure 1b). In the sMD simulations, the potential was scaled down uniformly, which makes the sMD simulations become equivalent to unbiased MD simulations of the same system at a high temperature. We therefore believe that the motion of loop C in the sMD simulations most likely reflects the natural rearrangement of the loop. In fact, loop C has been found to play important roles in acetylcholine receptor-channel gating and ligand binding 22. The functional study from Jadey et al. indicated that loop C of AChR has a “catch” and “hold” mechanism, in which it undergoes an “open” to a “closed” state conformational change during the binding process 23. From the alignment of the crystal structures of AChBP in the apo form and its complexes with agonists, partial agonists and antagonists, we also found that loop C is able to adopt different conformations (Figure 1d). We further measured the RMSDs of loop C using the sMD simulation data with respect to the crystal structures by using the trajectory from an sMD simulation (Figure S4). The result shows that before the ligand left the binding pocket, the conformation of loop C is closer to the “close” structure (PDB: 2BR7), while after the unbinding event, loop C adopts the conformations closer to the “open” structure (PDB: 2C9T). Thus, the conformational change of loop C observed in sMD simulations is validated by experiment, which further motivated us to consider this conformational change as an additional factor in the unbinding studies using InMetaD simulations.

Metadynamics with loop C “opening” and “closing” as an additional CV To investigate the role of loop C in the unbinding process, the progression of “opening” and “closing” of loop C as an additional path-CV (CVloop) was added to the previous CV set that involves only CVdis and CVsolv. Using this new set of CVs, a second round of twenty independent InMetaD simulations were carried out. The residence times thus obtained were then subject to the KS-test. We found that the p-value boosted from 0.003 to 0.15 (Figure 2a), which is about 50-fold improvement compared to that from the previous InMetaD simulations and reflects that the residence times follow satisfactorily a time-homogeneous Poisson distribution 13,14. Thus, the


Page 6 of 19

Page 7 of 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


second round of InMetaD simulations verified that the progression of “opening” and “closing” of loop C is indeed a critical CV closely related to the ASEM unbinding process. The mean residence time τ was obtained by fitting the ECDF to TCDF (see eq (5)) and the corresponding unbinding rate koff as the reciprocal of τ was calculated as 0.0033 s-1(see Method section), which is comparable to the experimental value of 0.0003 s-1. 6 It is notable that the experimental data was obtained from the in vivo PET tracer assay, which was measured in a more complicated biological environment rather than a simple protein-ligand system, and could be longer than the residence time of the ligand in the receptor. The low off-rate is also reflected by the high binding affinity of ASEM (~0.4 nM) 18, as a high binding affinity results in a stronger binding of the ligand to the receptor and makes it more difficult to dissociate 2. These two aspects verify indirectly the validity of our InMetaD simulations and that the selected CVs have captured the most relevant slow degrees of freedom for the ASEM unbinding process. Stable states during the unbinding process One of the important aims of molecular simulations is to reveal the mechanism of a microscopic process at the atomistic scale which can hardly be obtained from experiment, in the present case for the ASEM unbinding process. Firstly, we want to find the possible intermediates during the unbinding process. From the sMD and InMetaD simulations, two dominant states were identified as shown in Figure 2. These two states - state U and state D - can be characterized by the orientation of the sulfone group on ASEM: In state U the sulfone group points up; in state D, the sulfone group points down. The two states can remain stable for at least 200 ns simulations, as indicated by the low fluctuations of the RMSDs in Figure S1 and S2.

Figure 2. Cumulative distribution of the residence times from the second round of InMetaD



simulations, using the distance and solvation CVs, together with the loop C path CV, and the representative structures taken out from the simulations. (a) cumulative distribution of the residence times from 20 independent InMetaD simulations. (b-d) the representative structure of the bound state, state U and state D. ASEM and key residues are shown in stick form. Protein and loop C are shown in cartoon. Hydrogen bonds are shown as dashed lines. To estimate the residence times of the two states, twenty additional InMetaD simulations were conducted for each state using the same CV set and stopped once the system escaped from its initial state. The residence time was calculated by eq (4). As a comparison, the residence time of the bound state was also calculated. For the latter state, the time was taken from the previously full unbinding simulations but only counted the time until the system left the bound state instead of the time until it became fully unbound. Therefore, the calculated time corresponds to the residence time for the bound state rather than to the full residence time. The reweighted residence times and the p-value for each state from the KS-test are shown in Table S1. Of all the states identified, the bound state has a mean residence time around 180 s. It is far more stable than the two other states we identified. The calculated mean residence times of state U (138 ms) and state D (98 ms) are thus much shorter than that of the bound state. This result matches our previous study of ASEM in the sense that it showed that this bound state is a global minimum in the binding pocket 24. The residence times imply that the U and D states are much less stable than the bound state but they still show considerable stability. Thus, from the structural point of view, it is desirable to understand the molecular determinants of the stability of the two states. As described in the previous binding mode study 24, the protonated tip nitrogen of the diazobicyclic head group of ASEM forms a hydrogen bond with the backbone oxygen of Trp145 in the bound state (Figure 2). In the U and D states, the diazobicyclic group also stays inward and anchors ASEM in the binding site. In the bound state, the dibenzo[b,d]thiophene ring points to loop C and becomes π-π stacking with Trp53. However, neither state U nor state D can form such an interaction - in both states, Trp53 flips up and likely forms an edge-to-face stacking with the dibenzo[b,d]thiophene ring of ASEM. Previous study indicated that the π-π stacking interaction between Trp53 and dibenzo[b,d]thiophene ring makes an important contribution to the stabilization of the ligand. Compared to the edge-to-face mode, the parallel face-to-face π-π stacking gained ~4 kcal/mol more favorable energy 24. In state U and state D the ligand loses some favorable interactions with Trp53 due to the edge-to-face mode, but newly formed hydrogen bond interactions are able to stabilize ASEM in these two states. In state U, the sulfone group of ASEM rotates upward and forms hydrogen bonds with Gln55 and Gln114, while in state D, this group rotates downward and forms a hydrogen bond with Arg182. These hydrogen bonds stabilize ASEM and can explain why state U and state D can be identified from both InMetaD and sMD simulations and show stability during the conventional MD simulations. Unbinding mechanism and rate limiting steps The individual kinetic constants kij (Table S2), i.e. the transition rate from state i to state j, are calculated from the additional InMetaD simulations for each transition as described previously. From the state-to-state transition model, a Markov state model was built. As summarized in Figure 3, the ASEM unbinding process comprises of two steps. The first step is from the bound state to


Page 8 of 19

Page 9 of 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


state U or D, and the second step is from state U or D to the unbound state. Reverse transitions and the transitions between state U and state D were also found.

Figure 3. The state-to-state transition network, representative structures of the stable states, and representative structures for the transitions between the stable states. The bold rounded rectangle panels represent the stable states, i.e. the bound state, state U, and state D. The number below each panel is the residence time of the state. The pink rectangle panels are the representative intermediate structures between the stable states. B-U, B-D, U-D, U-Ub, and D-Ub indicate the representative intermediate structures between the bound state and state U (B-U), the bound state and state D (BD), state U and state D (U-D), state U and the unbound state (U-Ub), and state D and the unbound state (D-Ub), respectively. The arrow between two stable states represents the transition direction and the number next to the arrow is the transition rate (s-1). ASEM (cyan) and key residues (Trp53, shown in green) are shown in stick form. Protein (white) and loop C (magenta) are shown in cartoon. To evolve from the bound state to the unbound state, the system first needs to reach either state U or state D. The rate for such a transition is rather low – it was predicted to be 0.0885 s-1 from the bound state to state U and 0.0221 s-1 to state D. The rate from the bound state to state U is slightly higher than that to state D, indicating that the conformational changes from the bound state to state U is relatively easy. In both state U and state D, Trp53 flips up as compared to its conformation in the bound state. This conformational change would involve high energy conformations such as BU and B-D, which are representative structures extracted from the transitional trajectories between the bound state and state U (B-U) and those between the bound state and state D (B-D). The energy barrier from the bound state to state U is lower when the sulfone group on ASEM rotates up. This



conformational change has less steric conflicts with Tyr184 and Glu185 than that encountered during the transition from the bound state to state D. Once the system reaches state U or state D, ASEM becomes comparatively easier to dissociate. The rates from state U and state D to the unbound state are 7.2 s-1 and 64.0 s-1, respectively. While state D has a lower rate to reach the bound state, its evolution to the unbound state is easier (with a higher rate) than that from state U. The representative structures U-Ub and D-Ub were extracted from the unbinding trajectories starting from state U and from state D, respectively. As shown in Figure 2c and Figure 3, to dissociate from state U, ASEM needs to break the hydrogen bonds with Gln55 and Gln114 and cross the barrier formed by Glu185, Cys186 and Cys187. The conformational barrier ASEM needs to cross from state D has less steric conflicts, since in state U the sulfone group of ASEM has already pointed out to form a hydrogen bond with Arg182 that originally points out to the solvent. Once the hydrogen bond breaks, it flips back and ASEM leaves the binding site. Reversible transitions were also found between the states - state U and state D not only can evolve to the unbound state but also can return to the bound state. Compared to the transition rates from the bound state, the rates for the reverse processes are much higher - the rates from the two states to the bound state are 137.7 s-1 and 149.3 s-1, respectively. Furthermore, direct transitions between state U and state D were also found. The conformational change between state U and state D can occur without going through the bound state. Compared with the hydrogen bonds with Gln55 and Gln114 in state U, the hydrogen bond with Arg182 in state D is easier to break. The transition rate from state D to U (712.4 s-1) is approximately four-fold higher than that from state U to D (163.0 s-1). The conformational changes between U and D do not incur the high energy conformational flip of Trp53 and ASEM has larger space to rotate within the pocket when Trp53 is vertical, so the transitions between the two states can be frequently found. This also implies that the unbinding process may not proceed via simple routes such as from the bound state to state U, and then from state U to unbound. Once the system evolves to state U or state D there are more alternative paths for ASEM to unbind such as the bound state - state U - state D - unbound or bound - state D - state U – the unbound state. By solving the eigenvalues of the master equation associated with the kinetics matrix, slow steps in the unbinding network were identified. The slowest rate 0.012 s-1, which was obtained from the first non-zero eigenvalue, corresponds to the transition from the bound to the unbound state. This rate is in fair agreement with the overall rate calculated from the full unbinding metadynamics (0.0033 s1). The second eigenvalue and its associated vector (Figure S3) show that the transitions associated with state U constitute rate-limiting steps. Although the transition rate from the bound state to state U is higher than that to state D (0.0885 s-1 vs. 0.0221 s-1), the transition from state U to the unbound state is slower than that from state D to the unbound state (7.2 s-1 vs. 64.0 s-1). In addition, the eigenvector (Figure S3) shows that the most probable transition from state U is to the bound state, then to state D, and finally to the unbound state. From state U the system preferably returns to the bound state instead of evolving to the unbound state (137.7 s-1 vs. 7.2 s-1). Therefore, besides the bound state, state U constitutes a second minimum that traps the ligand and lowers the unbinding


Page 10 of 19

Page 11 of 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


rate. Discussion From the predicted unbinding mechanism of ASEM from the α7-AChBP chimera structure, we learn that ASEM leaves the binding pocket mainly via two steps. In the first step, the dibenzo[b,d]thiophene ring rotates to either the U state or the D state. Meanwhile Trp53 rotates from being vertical to horizontal. In this step loop C opens up a little bit to create sufficient space for the conformational changes from the bound state to the intermediate states. In the second step, ASEM is stabilized by forming the hydrogen-bond interactions with Gln55 and Gln114 in state U and with Arg182 in state D. However, since in state U and state D the interaction between the dibenzo[b,d]thiophene ring of ASEM and the indole ring of Trp53 is weaker than when Trp53 is in the vertical form in the bound state, ASEM can leave the binding pocket with a higher rate. Experimental evidence indicates that the binding of a ligand to the nAChR pocket may occur through an induced fit or a “catch-and-hold” mechanism 23. Our simulations of the unbinding process describe such an induced fit mechanism at the atomistic level, assuming that the unbinding process is the reverse of the binding process. In the reverse mechanism, ASEM probably encounters the α7 structure to form the U and D states with Trp53 in the vertical form. Then via conformational rearrangement ASEM forms a more stable bound state with Trp53 horizontally. In the binding/unbinding process, loop C plays an important role not only during the encountering events to “catch” the ligand but also during the “hold” process when it opens up a little bit to permit the conformational change from state U or D to the final bound state. The steps involved in forming the U or D state include: (1) the rotations of Trp53 and the dibenzo[b,d]thiophene ring of ASEM in the pocket during the transition from the bound state to state U or D; (2) formation of the hydrogen bonds between ASEM and Gln55 and Gln114, and; (3) the open up of loop C during the transition from the bound state to state U or D. In steps (1) and (2), the dibenzo[b,d]thiophene ring has to flip from the horizontal to the vertical form and meanwhile Trp53 rotates from the horizontal to the vertical form. The limited space in the binding pocket implies a high energy barrier and a low transition rate to induce such conformational changes. From the metadynamics trajectories obtained using only two CVs, we found that there are essentially no such conformational changes of the dibenzo[b,d]thiophene ring or the indole ring of Trp53 because loop C did not open to generate more space. Only in the trajectories that include loop C conformational changes, the opening of loop C can offer more space for such conformational changes. This explains why the conformational change of loop C is essential for the unbinding. In this study we used the α7-AChBP chimera structure instead of the full-length α7-nAChR membrane protein in order to simplify the model and also because of the mere fact that the fulllength α7-nAChR structure is not available yet. However, since the residues in the binding pocket of the α7-AChBP chimera structure are lined entirely by the α7-nAChR residues, the mechanisms and metastable states occurring in the AChBP most probably also exist in α7-nAChR and can therefore provide valuable information for the rational development of α7-nAChR PET tracers. The sulfone group of ASEM can in the binding pocket form hydrogen-bond interactions with Gln55 and Gln114 in state U and with Arg182 in state D. However, other α7-nAChR ligands do not have such



kind of sulfone group and will thus not be stabilized by this type of interactions. This explains why ASEM sticks out with excellent properties among the PET radioligands, something that can hardly be explained by considering solely the bound state. Selection of CVs in metadynamics is an essential step that often requires time consuming tries and tests. For a complex system in which CVs are difficult to identify, it is helpful to use sMD to sample its conformational rearrangement, which normally occurs at much longer timescales, in order to identify the key conformational changes so as to inspire the selection of CVs. In this work, we benefit from the application of sMD in identifying important CVs. However, our selection of the additional CV was mainly based on the inspection of the sMD trajectories. There are several methods developed for automatic selection of CVs, such as SGOOP and TICA 12,25–27. These methods can automatically adjust the reaction coordinate by a linear combination of a bunch of candidiate CVs. We believe that the combination of sMD with these methods could further facilitate the selection of CVs for InMetaD simulations. Conclusion In this work, we combined sMD and InMetaD simulation techniques to reveal the unbinding mechanism of ASEM from the α7-AChBP. In our initial InMetaD simulations of the unbinding kinetics of ASEM from the α7-AChBP chimera structure with the commonly used CVs (i.e. the distance CV and solvation CV), we have shown that it is not possible to obtain a statistically reasonable residence time. By applying sMD simulations, we identified that motion of loop C is critical for the ASEM unbinding process. When the progression of “opening” and “closing” of loop C was selected as an additional CV, the InMetaD simulation result indicates that the new CV set can reflect the most relevant slow degrees of freedom critical to the ASEM unbinding process. We found that there are more than one metastable state in the ASEM unbinding process and the identified rate limiting steps are highly associated with the interactions between the residues of the binding pocket and the sulfone group of ASEM. We also calculated the state-to-state transition rates and identified rate-limiting structures to provide a detailed unbinding mechanism of ASEM from the α7-AChBP chimera structure. sMD is CV free and can sample the conformational rearrangement of a protein-ligand system occurring at the timescale beyond seconds. This will help us identify CVs connected to the protein conformational changes critical to the ligand unbinding process. By carrying out InMetaD simulations with the sampling on such CVs also enhanced, we are able to reveal the mechanism of the unbinding process. Our strategy of combining sMD and InMetaD simulation techniques thus provides an effective way for studying the unbinding kinetics of protein-ligand systems in general.

Method Simulation setup The model of the α7-AChBP chimera structure in complex with ASEM was obtained from the previous study 24. The AMBER ff99SB-ILDN force field 28 and the general AMBER force field (GAFF) 29 were used for the protein and ASEM, respectively. The restrained electrostatic potential


Page 12 of 19

Page 13 of 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


(RESP) derived charges were used for ASEM with the electrostatic potential calculated at the Hartree-Fock level with the 6-31G* basis set using Gaussian 09 30. The TIP3P 31 water model was used to solvate the complex and 140 Na+ and 138 Cl- ions were used to neutralize the system. The Gromacs program 32 was employed for the MD simulations. Before production runs, energy minimization and a 200 ps restrained simulation in the NVT ensemble (T = 300 K) were carried out, followed by a 500 ps simulation carried out in the NPT ensemble (T = 300 K, P = 1 atm) for equilibrium. A time step of 2 fs was used in all the simulations and the LINCS algorithm 33 was applied to constrain the bonds involving hydrogen atoms. The Particle Mesh Ewald (PME) method 34 was used to recover the long range electrostatic interactions. The cut-offs for the short-range electrostatic interactions and van der Waals interaction were set to 10 Å. Infrequent metadynamics simulations Infrequent metadynamics (InMetaD) 20,35 was employed in this study. In an InMetaD simulation, a history-dependent biased potential V(s,t) is applied over a selected set of collective variables (CVs), with 2 ′ 𝑑 (𝑠𝑖 ― 𝑠𝑖(𝑡 )) 𝑉(𝑠,𝑡) = ∑𝑡′ < 𝑡ωexp ( ― ∑𝑖 = 1 2𝜎2 ), (1) 𝑖

where ω is the height and 𝜎𝑖 is the width of the bias potential deposited over the ith CV. In

this work, well-tempered method 36 was used and the Gaussian height was decreased according to ―𝑉(𝑠,𝑡′)

ω = ω0exp((1 ― 𝛾)𝑘𝐵𝑇),

(2)

where ω0 is the initial height, and 𝛾 is the bias factor that smoothly tunes the potential. The distance between the ligand and the binding pocket of the protein and the number of water molecules around the ligand have been considered as two CVs relevant to a ligand unbinding process and have been successfully used in several studies 13,14. In our initial InMetaD simulations, CVdis (i.e. the distance CV) was defined as the center of mass distance between the ligand (represented by the heavy atoms of ASEM) and the binding pocket (consisting of the Cα atoms of Ser32, Leu33, Ser34, Trp53, Leu54, Leu106, Leu116, Pro117, Ser144, Trp145, Ser146, Asp160, Tyr184, Glu185, Cys186, Cys187 and Tyr191), and CVsolv (i.e. the solvation CV) was defined as the number of water molecules (represented by the water oxygens) within 8 Å from the ligand center of mass. As CVs need to be smoothly differentiable 14, a switch function was used for CVsolv as described in the Supporting Information (eq. (S1)). A path CV 37–39 of loop C was used as an additional CV in our second round of InMetaD simulations. In a path CV, the progression along a reference path (Spath) reads 𝑛

𝐶𝑉𝑝𝑎𝑡ℎ =

∑𝑖 = 1𝑖𝑒𝑥𝑝( ― 𝜆𝑚𝑠𝑑(𝑹, 𝑹𝒊)) 𝑛

∑𝑖 = 1𝑒𝑥𝑝( ― 𝜆𝑚𝑠𝑑(𝑹, 𝑹𝒊))

,

(3)

where R denotes the coordinates during a simulation, Ri are the coordinates that compose the path, msd(R, Ri) is the mean-square deviation of R and Ri , and λ is a smoothing parameter the value of which is roughly proportional to the inverse of the average mean square displacement between two successive frames along the reference path. In our study, the Spath for the conformational change of



Page 14 of 19

loop C was defined from the representative crystal structures of AChBP with different loop C conformations (PDB: 2BR7, 2C9T, 2Y54, 2XYT, and 5XGL). 5 frames with uniform nearest neighbor fluctuation (i.e. mean-square deviation difference) were obtained using Morph Server and were used to describe the loop C conformational change. Note we only used the Cα atoms from Ala177 to Val200 as the reference and no bias was deposited in the orthogonal directions, since the distance from the Spath , which is defined as Zpath, is limited and the unbiased thermo-fluctuation could cross the energy barrier in this direction. The time evolutions of Spath and Zpath in the second round of twenty independent InMetaD simulations are shown in Figure S5. Plumed 2.3 40,41 patched Gromacs was used for the InMetaD simulations. The initial Gaussian height was 0.2 kcal/mol. One of the key assumptions for this method is that no bias is added in the transition state region. To achieve this, the frequency of the bias deposition should be slow. The pace for the bias deposition was set to 10 ps, which is the same as that used in previously InMetaD studies 13,14,20. From the InMetaD simulation time tMTD, the residence time t can be reweighted as: t = t𝑀𝑇𝐷𝛼 , (4)

〈𝑒

where α is the acceleration factor

𝑉(𝑠, 𝑡𝑀𝑇𝐷) k 𝐵T

〉, and s represents the CVs. We performed twenty

independent InMetaD simulations. The residence times obtained according to eq.(4) were used to build an empirical cumulative distribution function (ECDF) which was then fit to a theoretical cumulative distribution function (TCDF) to obtain the characteristic residence time τ: 𝑡

TCDF = 1 ― e

―(𝜏)

(5)

Finally, the difference between ECDF and the TCDF were assessed by the two sample Kolmogorov−Smirnov (KS) test. The null hypothesis was accepted if the p-value from the KS test was higher than the significance threshold (0.05). Potential-scaled MD simulations In an sMD simulations, the potential constructed from a force field is uniformly scaled by a factor β that is smaller than one. After tuning different β values and referring to previously studies 15,16, we scaled the potential energy by a factor of 0.4 (i.e. β=0.4) to compromise a reasonable simulation time. Heavy atoms of the protein backbone were restrained with a weak harmonic potential (with the force constant k = 50 kJ mol-1 nm-2) except for residues within 6Å of the ligand ASEM. In this way, the overall protein structure is preserved while the rearrangement of the binding site is allowed to be sampled during the ligand unbinding process. The unrestrained residues are: Ser32-Ser34, Leu36, Phe52-Gln55, Ala89, Tyr91, Thr101, Pro102, Leu106, Leu116-S118, Gln143-His148, Glu158-Asp160, Ser162, Gly163, Arg182-Asp193, and Phe196. Twenty simulation runs were performed using the Gromacs program with the scaled force field. Each simulation was stopped once the ligand was fully unbound from the pocket. State identification and construction of the Markov Model To provide details of the unbinding mechanism, a Markov Model as described by Tiwary et al. was built 13,14. Intermediate states where the system spends much longer time than average during the unbinding process were identified by conformational clustering. To maximally avoid possible artifacts generated from the sMD or InMetaD simulations, we only considered the states that can be


Page 15 of 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


found from both types of simulations and neglected those which can only be identified from one type of simulations. Firstly, protein heavy atoms were aligned with the trajectories obtained from the sMD and InMetaD simulations with the bias added over CVdis, CVsolv, and CVloop. Next, the states in which the system spends a minimum of 5 ns in the InMetaD simulations and 1 ns in the sMD simulations were considered and a cluster analysis based on the ligand RMSDs was performed with a cutoff of 0.2 nm. Conformations were taken out from the clusters most populated in sMD and InMetaD simulations and were further verified by 200 ns unbiased MD simulations. Although this procedure may ignore some intermediate states, as discussed in the result section, the most populated intermediates can be identified with minimum artifact. For each of the states verified via the unbiased MD simulations, twenty independent InMetaD simulations with the bias added over CVdis, CVsolv, and CVloop were performed to obtain the residence time. Each simulation was stopped once the system escapes from the starting state. The KS test was performed for each set of simulations. For each transition, the rate was calculated as described in refs [10] and [11]. The transition kinetic matrix was built and the rate-determining step was obtained from the most negative eigenvalue of the transpose matrix. Supporting Information Supporting Information Available: The definition of the switch function for CVsolv. Figure S1: Evolution of RMSD for ASEM in state U. Figure S2: Evolution of RMSD for ASEM in state D. Figure S3: Eigenvector of the kinetic transition matrix. Figure S4: RMSDs of loop C along an sMD simulation with respect to the crystal structures. Figure S5: Time evolution of Spath and Zpath along the second round of the twenty independent InMetaD simulations. Table S1: Residence time and pvalue for each state obtained from the KS-test. Table S2: Kinetics transition matrix. (PDF) Acknowledgment The authors acknowledge support from the Swedish Foundation for Strategic Research (SSF) through the project “New imaging biomarkers in early diagnosis and treatment of Alzheimer's disease”. We would like to thank Dr. Yong Wang for his helpful discussions and comments. Y.Z acknowledges financial support from the China Scholarship Council (CSC) for his PhD studies. The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC.

Competing interests The authors declare that there are no competing interests.

References (1) Tonge, P. J. Drug-Target Kinetics in Drug Discovery. ACS Chem. Neurosci. 2018, 9, 29–39. https://doi.org/10.1021/acschemneuro.7b00185. (2) Copeland, R. A. The Drug-Target Residence Time Model: A 10-Year Retrospective. Nat. Rev. Drug Discovery. 2016, 15, 87–95. https://doi.org/10.1038/nrd.2015.18. (3) Bernetti, M.; Cavalli, A.; Mollica, L. Protein–Ligand (Un)Binding Kinetics as a New Paradigm



for Drug Discovery at the Crossroad between Experiments and Modelling. Med. Chem. Commun. 2017, 8, 534–550. https://doi.org/10.1039/C6MD00581K. (4)

Zhang, R.; Monsma, F. Binding Kinetics and Mechanism of Action: Toward the Discovery and Development of Better and Best in Class Drugs. Expert Opin. Drug Discovery. 2010, 5, 1023– 1029. https://doi.org/10.1517/17460441.2010.520700.

(5)

de Witte, W. E. A.; Danhof, M.; van der Graaf, P. H.; de Lange, E. C. M. In Vivo Target Residence Time and Kinetic Selectivity: The Association Rate Constant as Determinant. Trends Pharmacol. Sci. 2016, 37, 831–842. https://doi.org/10.1016/j.tips.2016.06.008.

(6)

Hillmer, A. T.; Li, S.; Zheng, M. Q.; Scheunemann, M.; Lin, S. fei; Nabulsi, N.; Holden, D.; Pracitto, R.; Labaree, D.; Ropchan, J.; Teodoro, R.; Deuther-Conrad, W.; Esterlis, I.; Cosgrove, K; Brust, P.; Carson, R.; Huang, Y. PET Imaging of Α7 Nicotinic Acetylcholine Receptors: A Comparative Study of [18F]ASEM and [18F]DBT-10 in Nonhuman Primates, and Further Evaluation of [18F]ASEM in Humans. Eur. J. Nucl. Med. Mol. Imaging. 2017, 44, 1042–1050. https://doi.org/10.1007/s00259-017-3621-8.

(7)

Votapka, L. W.; Amaro, R. E. Multiscale Estimation of Binding Kinetics Using Brownian Dynamics, Molecular Dynamics and Milestoning. PLoS Comput. Biol. 2015, 11, e1004381. https://doi.org/10.1371/journal.pcbi.1004381.

(8)

Husic, B. E.; Pande, V. S. Markov State Models: From an Art to a Science. J. Am. Chem. Soc. 2018, 140, 2386–2396. https://doi.org/10.1021/jacs.7b12191.

(9)

Zuckerman, D. M.; Chong, L. T. Weighted Ensemble Simulation: Review of Methodology, Applications, and Software. Annu. Rev. Biophys. 2017, 46, 43–57. https://doi.org/10.1146/annurev-biophys-070816-033834.

(10)

Pang, X.; Zhou, H.-X. Rate Constants and Mechanisms of Protein–Ligand Binding. Annu. Rev. Biophys. 2017, 46, 105–130. https://doi.org/10.1146/annurev-biophys-070816-033639.

(11)

Tiwary, P.; Limongelli, V.; Salvalaglio, M.; Parrinello, M. Kinetics of Protein–Ligand Unbinding: Predicting Pathways, Rates, and Rate-Limiting Steps. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, E386–E391. https://doi.org/10.1073/pnas.1424461112.

(12)

Tiwary, P. Molecular Determinants and Bottlenecks in the Dissociation Dynamics of BiotinStreptavidin. J. Phys. Chem. B. 2017, 121, 10841–10849. https://doi.org/10.1021/acs.jpcb.7b09510.

(13)

Casasnovas, R.; Limongelli, V.; Tiwary, P.; Carloni, P.; Parrinello, M. Unbinding Kinetics of a P38 MAP Kinase Type II Inhibitor from Metadynamics Simulations. J. Am. Chem. Soc. 2017, 139, 4780–4788. https://doi.org/10.1021/jacs.6b12950.

(14)

Tiwary, P.; Mondal, J.; Berne, B. J. How and When Does an Anticancer Drug Leave Its Binding Site? Sci. Adv. 2017, 3, e1700014 31. https://doi.org/10.1126/sciadv.1700014.

(15)

Mollica, L.; Decherchi, S.; Zia, S. R.; Gaspari, R.; Cavalli, A.; Rocchia, W. Kinetics of ProteinLigand Unbinding via Smoothed Potential Molecular Dynamics Simulations. Sci. Rep. 2015, 5, 11539. https://doi.org/10.1038/srep11539.

(16)

Mollica, L.; Theret, I.; Antoine, M.; Perron-Sierra, F.; Charton, Y.; Fourquez, J. M.; Wierzbicki, M.; Boutin, J. A.; Ferry, G.; Decherchi, S.; Bottegoni, G.; Ducrot P.; Cavalli, A. Molecular Dynamics Simulations and Kinetic Measurements to Estimate and Predict ProteinLigand Residence Times. J. Med. Chem. 2016, 59, 7167–7176. https://doi.org/10.1021/acs.jmedchem.6b00632.

(17)

Schuetz, D. A.; Bernetti, M.; Bertazzo, M.; Musil, D.; Eggenweiler, H. M.; Recanatini, M.;


Page 16 of 19

Page 17 of 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


Masetti, M.; Ecker, G. F.; Cavalli, A. Predicting Residence Time and Drug Unbinding Pathway through Scaled Molecular Dynamics. J. Chem. Inf. Model. 2019, 59, 535–549. https://doi.org/10.1021/acs.jcim.8b00614. (18)

Horti, A. G. Development of [18 F]ASEM, a Specific Radiotracer for Quantification of the A7NAChR with Positron-Emission Tomography. Biochem. Pharmacol. 2015, 97, 566–575. https://doi.org/10.1016/j.bcp.2015.07.030.

(19)

Callegari, D.; Lodola, A.; Pala, D.; Rivara, S.; Mor, M.; Rizzi, A.; Capelli, A. M. Metadynamics Simulations Distinguish Short- and Long-Residence-Time Inhibitors of CyclinDependent Kinase 8. J. Chem. Inf. Model. 2017, 57, 159–169. https://doi.org/10.1021/acs.jcim.6b00679.

(20)

Tiwary, P.; Parrinello, M. From Metadynamics to Dynamics. Phys. Rev. Lett. 2013, 111, 1–5. https://doi.org/10.1103/PhysRevLett.111.230602.

(21)

Tiwary, P.; Mondal, J.; Morrone, J. A.; Berne, B. J. Role of Water and Steric Constraints in the Kinetics of Cavity–Ligand Unbinding. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 12015–12019. https://doi.org/10.1073/pnas.1516652112.

(22)

Purohit, P.; Auerbach, A. Loop C and the Mechanism of Acetylcholine Receptor-Channel Gating. J. Gen. Physiol. 2013, 141, 467–478. https://doi.org/10.1085/jgp.201210946.

(23)

Jadey, S.; Auerbach, A. An Integrated Catch-and-Hold Mechanism Activates Nicotinic Acetylcholine Receptors. J. Gen. Physiol. 2012, 140, 17–28. https://doi.org/10.1085/jgp.201210801.

(24)

Kuang, G.; Zhou, Y.; Zou, R.; Halldin, C.; Nordberg, A.; Långström, B.; Ågren, H.; Tu, Y. Characterization of the Binding Mode of the PET Tracer [18F]ASEM to a Chimera Structure of the Α7 Nicotinic Acetylcholine Receptor. RSC Adv. 2017, 7, 19787–19793. https://doi.org/10.1039/C7RA00496F.

(25)

Tiwary, P.; Berne, B. J. Spectral Gap Optimization of Order Parameters for Sampling Complex Molecular Systems. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, 2839–2844. https://doi.org/10.1073/pnas.1600917113.

(26)

Pramanik, D.; Smith, Z.; Kells, A.; Tiwary, P. Can One Trust Kinetic and Thermodynamic Observables from Biased Metadynamics Simulations?: Detailed Quantitative Benchmarks on Millimolar Drug Fragment Dissociation. J. Phys. Chem. B. 2019, 123, 3672–3678. https://doi.org/10.1021/acs.jpcb.9b01813.

(27)

M. Sultan, M.; Pande, V. S.; Sultan, M. M.; Pande, V. S.; M. Sultan, M.; Pande, V. S. TICAMetadynamics: Accelerating Metadynamics by Using Kinetically Selected Collective Variables. J. Chem. Theory Comput. 2017, 13, 2440–2447. https://doi.org/10.1021/acs.jctc.7b00182.

(28)

Wickstrom, L.; Okur, A.; Simmerling, C. Evaluating the Performance of the FF99SB Force Field Based on NMR Scalar Coupling Data. Biophys. J. 2009, 97, 853–856. https://doi.org/10.1016/j.bpj.2009.04.063.

(29)

Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157–1174. https://doi.org/10.1002/jcc.20035.

(30)

M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ort, and D. J. F.



Gaussian 09. Gaussian, Inc.: Wallingford CT 2016. (31)

Price, D. J.; Brooks, C. L. A Modified TIP3P Water Potential for Simulation with Ewald Summation. J. Chem. Phys. 2004, 121, 10096–10103. https://doi.org/10.1063/1.1808117.

(32)

Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindah, E. Gromacs: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers. SoftwareX. 2015, 1–2, 19–25. https://doi.org/10.1016/j.softx.2015.06.001.

(33)

Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463–1472. https://doi.org/10.1002/(SICI)1096-987X(199709)18:123.0.CO;2-H.

(34)

Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577–8593. https://doi.org/10.1063/1.470117.

(35)

Wang, Y.; Valsson, O.; Tiwary, P.; Parrinello, M.; Lindorff-Larsen, K. Frequency Adaptive Metadynamics for the Calculation of Rare-Event Kinetics. J. Chem. Phys. 2018, 149, 072309. https://doi.org/10.1063/1.5024679.

(36)

Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 020603. https://doi.org/10.1103/PhysRevLett.100.020603.

(37)

Branduardi, D.; Gervasio, F. L.; Parrinello, M. From A to B in Free Energy Space. J. Chem. Phys. 2007, 126, 054103. https://doi.org/10.1063/1.2432340.

(38)

Díaz Leines, G.; Ensing, B. Path Finding on High-Dimensional Free Energy Landscapes. Phys. Rev. Lett. 2012, 109, 020601. https://doi.org/10.1103/PhysRevLett.109.020601.

(39)

Zhou, Y.; Hussain, M.; Kuang, G.; Zhang, J.; Tu, Y. Mechanistic Insights into Peptide and Ligand Binding of the ATAD2-Bromodomain via Atomistic Simulations Disclosing a Role of Induced Fit and Conformational Selection. Phys. Chem. Chem. Phys. 2018, 20, 23222–23232. https://doi.org/10.1039/c8cp03860k.

(40)

Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello M. PLUMED: A Portable Plugin for FreeEnergy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009, 180, 1961– 1972. https://doi.org/10.1016/J.CPC.2009.05.011.

(41)

Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New Feathers for an Old Bird. Comput. Phys. Commun. 2014, 185, 604–613. https://doi.org/10.1016/J.CPC.2013.09.018.


Page 18 of 19

Page 19 of 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


For Table of Contents Only


Enhanced sampling simulations of ligand unbinding kinetics

Recommend Documents