In Silico Reoptimization of Binding Affinity and ... - ACS Publications

Jul 16, 2019 - central aromatic ring of RL-45 displays immaculate flexibility in bypassing any steric interaction resulting due to T338M mutation. Acc...
0 downloads 0 Views 1MB Size
Subscriber access provided by KEAN UNIV

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

In Silico Re-Optimization of Binding Affinity and Drug-Resistance Circumvention Ability in Kinase Inhibitors:Case study with RL-45 and Src kinase Jaya Krishna Koneru, Suman Sinha, and Jagannath Mondal J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.9b02883 • Publication Date (Web): 16 Jul 2019 Downloaded from pubs.acs.org on July 17, 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 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

In silico Re-optimization of Binding Affinity and Drug-resistance Circumvention Ability in Kinase Inhibitors: Case study with RL-45 and Src kinase Jaya Krishna Koneru, Suman Sinha* and Jagannath Mondal* Tata Institute of Fundamental Research, Centre for Interdisciplinary Sciences, 36/P Gopanpally, Serilingampally Mandal, Hyderabad, 500107, India *Email: [email protected], [email protected]

Abstract: A major bottle neck in the development of kinase inhibitors has been the onset of drug resistance around the gate-keeper residues of Src kinase. Although recent time has seen the reports of certain second generation kinase inhibitors which are capable of bypassing the drug resistance by circumventing kinase mutation, their kinase binding efficacy has remained considerably weaker than the classical ATP-competitive kinase inhibitors. Using a recently synthesized second generation kinase inhibitor RL-45 as a template, the current work integrates fragment-based drug discovery and quantitative structure activity relationship study with enhanced molecular dynamics simulation approaches, namely metadynamics and replica exchange free energy perturbation, and demonstrates how one can optimally redesign and assess novel Src kinase inhibitors, by minimal introduction of new fuctional moeities around template kinase inhibitor. Interestingly, unlike many synthetic kinase inhibitors, these in-silico optimized small molecule derivatives of RL-45 are found to be potentially capable of serving dual purposes, crucial for efficacy of an ideal kinase inhibitor: namely, a) circumventing gate-keeper residue mutation related drug resistance in Src kinase, unlike many commercial kinase inhibitors and b) manifesting superior resilience against unbinding from the kinase active site. The computer simulation, boosted by enhanced sampling techniques, further reveals that these designed inhibitors bring about key interactions in the form of significantly long-standing hydrogen bonds and hydrophobic pocket, otherwise weak in the template bioactive kinase inhibitor, which enhance the binding efficacy of these newly designed ligands in the kinase binding pocket.

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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: Protein kinases are defined by their ability to catalyze the transfer of the terminal phosphate of ATP to substrates that usually contain a serine, threonine or tyrosine residue. Hence kinases are known for controlling the transfer of energy and regulating many aspects of cell life. They typically share a conserved arrangement of secondary structure elements, including a well-conserved activation loop, which is important in regulating kinase activity and is marked by conserved aspartic acid-phenylalanine-glycine (DFG) motif at the end of the loop.1,2 Abnormal kinase activity correlates with many diseases including cancer and therefore development of drugs, which regulate the aberrations in kinase activity (i.e. socalled kinase inhibitors), are in full swing3. Receptor tyrosine kinase signalling pathways have been successfully targeted to therapeutically intervene for cancer treatment with numerous successful examples. In many cases, kinase inhibitors compete with the ATP to bind to the designated cavity and essentially stall the activity of kinase.4 Although more than a dozen kinase inhibitors are Food and Drug Administration (FDA)-approved, with several more in clinical trials, the onset of drug resistance remains a fundamental challenge in the development of these kinase inhibitors.5 In a large fraction of cases, a major source of drug resistance is related to mutations in the targeted kinase6. The most common mutations occur at the gatekeeper residue in the hinge region of the kinase, and these directly prevent or weaken the interaction with the inhibitor7. This can be attributed to the fact that most kinase inhibitors are ATP-competitive molecules, called type I inhibitors4,6. In the most typical of these mutations, a relatively smaller amino acid side chain such as threonine (Thr) is exchanged for a larger hydrophobic residue such as methionine (Met) or isoleucine (Ile). Many type I kinase inhibitors including dasatinib have

2 ACS Paragon Plus Environment

Page 2 of 25

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

The Journal of Physical Chemistry

been found to be ineffective against the commonly found mutation of Thr to Met at the 338th residue in Src and Abl kinases (the so-called T338M mutation)7–10 Hence, a big thrust of the pharmaceutical industry has been the redesign of kinase inhibitors, which have the potential of circumventing the drug-resistance due to gate-keeper residue mutation. Accordingly, structure-based drug discovery strategies7,11–13 have recently been employed in designing type II hybrid kinase inhibitors which not only bind in the ATP pocket of kinases, but also extend past the gatekeeper residue into a less conserved adjacent allosteric site. The efforts have been successful in designing kinase inhibitors, which can efficiently circumvent resistance related to T338M gate-keeper residue mutation. Recent computational efforts by Mondal and coworkers14 have also provided a free-energetic basis and interpretation of a FDA approved drug dasatinib’s lost activity and retention of activity by one such designed inhibitor RL-457 in the face of adverse T338M gate-keeper mutation in src-kinase. However, a closer comparison of intrinsic kinase-binding affinity of such newly designed hybrid kinase inhibitor (which are able to avoid gate-keeper residue mutation) such as RL-45, with that of a FDA approved drug dasatinib (which is vulnerable to gate-keeper residue mutation) shows that the intrinsic binding potency of RL-45 to the wild-type kinase inhibitor remains several order of magnitude weaker (IC50=0.021 +/- 0.005 μM) relative to that of dasatinib (IC50=0.0004 +/- 0.0002 μM).7,14 In this context, the current work addresses a pertinent question: how can we optimize the kinase inhibitor’s ability of both withstanding the gate-keeper residue mutation as well as their stronger binding affinity towards kinase? Towards this end, we offer a robust computational approach (see Figure 1) to redesign the type II kinase inhibitor RL-45, such that the intrinsic binding affinity towards a wild-type srckinase is enhanced, while simultaneously retaining its ability to circumvent kinase gatekeeper residue mutation. Specifically, we integrate in-silico fragment-based drug discovery (FBDD) technique15, generally used to identify novel bioactive compounds with diverse 3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

scaffolds and quantitative structure-activity relationship (QSAR)16 with the enhanced sampling based free energy simulation technique, namely metadynamics17 and free energy perturbation/replica exchange with solute tempering (FEP/REST)18,19 to cross-validate the resilience of a number of small molecules against gate-keeper mutation-based drug resistance and unbinding from the kinase cavity.

Figure 1: The schematics of computational approach employed in the current work is represented. A. The layout of the flowchart highlighting various components of structure prediction processes for kinase redesign is shown. B. The QSAR correlation plot, incorporated in the ligand design process, with the statistics is highlighted. C. Also shown is the crystallographic pose of the parent ligand RL-45 inside the src-kinase. D. Finally, the two new redesigned kinase inhibitors are shown.

4 ACS Paragon Plus Environment

Page 4 of 25

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

The Journal of Physical Chemistry

Method and Materials: The method follows an integrative computational approach as present in Figure 1(A-D). Here we present the details of each of the methods present in the flow chart. Fragment-based drug design: Fragment based drug design15 (FBDD) protocol has been employed for the purpose of designing new candidates for screening as kinase inhibitor. In this approach, firstly the core fragment structures are derived from a bioactive molecule. In the current case, earlier developed kinase inhibitor RL-45 (PDB: 3F3W)7 served as the parent bioactive molecule, which was deconstructed and redesigned to obtain new kinase inhibitor candidates. ACFIS webserver20 was used in this purpose, which performs the cycle of deconstruction of bioactive molecule into fragments, fragment reassembly, relinking the fragments to the junction of the core fragment, energy minimization of the newly reconstructed ligands, and finally re-ranking the hit candidates by short MD simulations followed by molecular-mechanics/poisson-boltzman with solvent accessible surface area (MM/PBSA) calculation. Our designed ligands have been grown fragment-wise to fit into the previously unexplored regions of RL-45 binding pocket through novel modifications. The molecular attachments were chosen in such a way that the core of RL-45 remains intact and only specific functional groups were added to gain additional binding of the target. A total of 200 molecules were generated which were further validated using both visual inspection and in house quantitative structure activity relationship (QSAR) analysis.16 QSAR:

All

the

two-dimensional

descriptors

(thermodynamic,

spatial,

electronic

and topological parameters) were calculated for QSAR analysis using Vlife MDS software.16 Thermodynamic parameters describe free energy change during drug receptor complex formation. Spatial parameters are the quantified steric features of drug molecules required for 5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

its complimentary fit with receptor and electronic parameters describe non-covalent bonding between drug molecules and receptor. The data set was divided into training and test set using sphere exclusion method. Training set (1020 molecules) and test set (347 molecules) were respectively used to develop the QSAR model and assess the predictive power of the model, which is not included in model generation. A k-nearest neighbour based QSAR model16 was generated using 1367 molecules compiled from ChEMBL21 with the statistics shown inside Figure 1B and 21 ligands were first sort-listed for further studies (see Figure S1 for the chemical structures). Only those compounds, which showed greater pIC50 than both RL-45 and dasatinib, were decided to be the ones for further investigation. By this approach, we could re-optimize/validate ACFIS results from the perspective of ligand descriptors. In the subsequent steps, resultant compounds were systematically visually inspected of binding poses, size of the compounds (we did not want to incorporate bigger molecules given the fact that it could lead to synthetic as well as ADME complexity) and ligand-protein interactions were further considered to zero down on two compounds. As detailed below, the complex soobtained was subjected to enhanced molecular dynamics simulations to verify the resilience of these two designed ligands against biased exit from the binding pocket of the src kinase and their efficacy in presence of adverse T338M mutation. The complete QSAR dataset (in sdf format) is provided in the supporting information in the form of a zip file. Screening based upon molecular dynamics simulations: In this work, we performed two different variants of enhanced sampling molecular dynamics (MD) simulation: well-tempered metadynamics22 simulation to assess the stability of the computationally designed proteinligand complexes of RL-45 and two designed ligand molecules and free energy perturbation/replica exchange with solute tempering (FEP/REST)18 to test their relative ability to circumvent the drug resistance related to gate-keeper residue mutations.

6 ACS Paragon Plus Environment

Page 6 of 25

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

The Journal of Physical Chemistry

Models and parameters for MD simulation: All molecular dynamics simulation was performed using Charmm36 force field available for protein23. However, for modeling the ligands, we found that the Charmm-Generalized forcefield (CGenFF)-derived parameters24, which is routinely used, predicted high charge and van der waals penalty for each of these molecules. Accordingly we went about optimizing the force field parameters for each of the three ligands. In this regard, we used “General Automated Atomic Model Parameterization (GAAMP)” algorithm proposed by Huang and co-workers25 to develop the optimized force fields for RL-45 and each of the two redesigned molecules. Briefly, using this approach, ab initio quantum mechanical data obtained using HF/6-31G* basis set were used as target data. Atomic partial charges were then optimized by fitting of target quantum mechanical electrostatic potential. The water-solute interactions with hydrogen-bonding donor and acceptor groups were also fitted with quantum-mechanical data to further optimize the parameters. Subsequently, dihedral angles with low energy barriers (so-called soft dihedrals) were fitted into quantum mechanically derived dihedral angle distributions. The protein-ligand complex is solvated and charge-neutralized in a water box with TIP3P26 water molecules. The system is then minimized to its lowest possible energy followed by NVT equilibration at constant room temperature (𝑇 = 298𝐾). Initially, wild type src kinase, individually in complex with RL-45 and two other designed ligands, is subjected to 100 ns of equilibrium MD simulation. During the period of the simulation, the ligands were found to reside in their designated binding pocket. We have found that these soft dihedrals would play an important role in retaining the flexibility in these drug-like molecules required to avoid repercussions of gate-keeper residue mutation. The final snapshot of the production run is then used as initial configuration for subsequent meta-dynamics simulations. Metadynamics simulation: In the current work, well-tempered metadynamics simulation22 was employed to test the resilience of the computationally designed protein-ligand complex 7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

constituting the src-kinase and each of the two re-designed ligands and parent kinase inhibitor, namely RL-45, LIG-1028 and LIG-1346. The metadynamics simulations are practical alternatives of long unbiased MD simulations and are designed to allow efficient back-and-forth movement across large free energy barriers, and therefore sampling of the relative stability of different poses with practical computing resources, while still maintaining full atomistic resolution. The key idea in the use of metadynamics is to build a timedependent bias as a function of the chosen collective variable that samples the ligand movements in and around its binding pose. As the correct pose has (to the extent that it reproduces key features of the native structure) a stronger binding affinity, a properly calibrated metadynamics bias should result in preferred displacement of incorrect poses, while leaving the correct pose relatively stable. In this work, we employ distance between centre of mass of binding pocket and ligand, torsion angle between ligand and a stable helix near binding pocket as the collective variables along which we build a time-dependent bias. The employment of dual collective variables of distance and angle in this work is mainly motivated by earlier works of Parrinello and coworkers27,28 in sampling free energy landscapes of various protein-ligand complexes. The end result of a typical metadynamics simulation is an estimate of the underlying Boltzmann probability distribution or equivalently the potential of mean force, also called free energy. The bias V(S, t) is typically constructed in the form of periodically added repulsive Gaussians, where S is the chosen CV which could be multidimensional. At any time t, the free energy F(S) can be obtained from the deposited bias V(S, t) as per the following equation:

𝑉(𝑆,𝑡→∞) = ―

𝛥𝑇 𝐹(𝑆) + 𝐶 𝑇 + 𝛥𝑇

Here T is the simulation temperature (300 K), ΔT is the tempering factor through which the amplitude of the bias deposited at a point s in the collective variable space is tuned down with 8 ACS Paragon Plus Environment

Page 8 of 25

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

The Journal of Physical Chemistry

the number of times the system visits a given point, and C(t) is a time-dependent constant which is irrelevant for the present work since we are interested only in normalized probability distributions resulting from the free energy. Other numerical parameters include the initial Gaussian hill height h=1.2 kJ/mol, bias deposition interval τ=400 steps (0.8 ps), and Gaussian width w=0.017 nm (distance), 0.35 degree (Torsion). All simulations were performed using Gromacs29 and its PLUMED30 plugin. Seminal works by Copeland and coworkers31 had earlier demonstrated drug residence times inside the protein active site as the true metric of drug efficacy. This hypothesis, together with the requirement of ligand’s thermodynamic stability at binding pocket, has worked as a basis of the current work. However, as known through multiple previous studies, converging the full free energy landscape, as well as computing ligand residence times in metadynamics simulations, while feasible32–34, would be computationally very intensive due to requirement of assessment via poisson statistics35 and will be redundant in absence of any experimental off-rate values and binding free energies (for newly redesigned ligands). Rather, similar to Clark et al17, we decide to estimate the ‘average metadynamics time’ it takes to unbind the ligand from the protein-ligand complex, as a practical affordable metric. Accordingly, we run three sets of independent metadynamics simulations for each protein-ligand complex, each starting from the same cavity. However, as will be presented in results and discussion section, we have also validated the robustness of our metric via computation of approximate free energy profiles and via variation of metadynamics bias deposition frequency. FEP/REST simulation: A crucial aim of the current work is to precisely assess the ability of each of the two proposed ligands to circumvent the T338M gate-keeper residue mutation in src-kinase. In this regard, FEP/REST simulation18 was performed to quantify the relative binding affinity of RL-45 and each of the two designed ligands to the wild type and T338M mutated src-kinase. We employ a similar approach as in Mondal et al’s previous work14. Two 9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

different systems were simulated: one where the protein bound to ligand undergoes the T338M mutation (‘protein-in-complex’ state and the free energy change referred to as 𝛥 𝐺𝑐𝑜𝑚𝑝𝑙𝑒𝑥) and other system where the ligand-free-protein undergoes the same mutation (‘protein-in-solvent’ state and the free energy change referred to as 𝛥𝐺𝑠𝑜𝑙𝑣𝑒𝑛𝑡). The relative binding free energy for the binding of the ligand to the wild type versus the mutant is then 𝛥𝛥𝐺. 𝛥𝛥𝐺 = 𝛥𝐺𝑐𝑜𝑚𝑝𝑙𝑒𝑥 ― 𝛥𝐺𝑝𝑟𝑜𝑡𝑒𝑖𝑛 The FEP simulation was carried out in discrete steps for a series of coupling parameters 𝜆 with values ranging from 𝜆 = 0 (corresponding to the protein with Thr338) to 𝜆=1 (corresponding to the the protein with Met338). In this work, we use the heavy atoms of the 338th residue of the protein (the location of mutation) to define the hot region described in the beginning of this section. The topology parameters defining the hot region were generated using PMX utility developed by de Groot and coworkers36. The crystallographic structure (pdb: 3F3W) of wild-type kinase was chosen as the configuration corresponding to λ =0. Subsequently, the configuration corresponding to each of the rest of the λs (including all intermediate λ and λ=1 i.e. the T338M mutant) was prepared, in a sequential manner, by starting with an equilibrated configuration corresponding to previous λ and equilibrating the configuration by performing equilibrium MD simulations using appropriately scaled Hamiltonian (using pmx utility) for the specific λ. We use a total of 24 𝜆 windows for the FEP/REST simulations with the effective temperature profiles 300K, 340K, 386K, 437K, 496K, 563K, 639K, 725K, 822K, 933K, 1058K, 1200K, 1200K, 1058K, 933K, 822K, 725K, 639K, 563K, 496K, 437K, 386K, 340K and 300K respectively, ensuring same (300K) temperature for 𝜆=0.0 and 𝜆=1.0. Each 𝜆 window was sampled for a period of 5 ns for both the complex and solvent simulations using NPT ensemble conditions.

10 ACS Paragon Plus Environment

Page 10 of 25

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

The Journal of Physical Chemistry

Results and Discussion: The second-generation hybrid kinase inhibitor RL-45 (PDB: 3F3W) synthesized by Getlik and coworkers7 acts as the starting point of our interest for our drug discovery protocol (See Figure 1). It was shown earlier that while FDA-approved drug dasatinib loses its kinasebinding efficiency in presence of T338M gate-keeper residue mutation, a new kinase inhibitor, coined RL-457, could circumvent the T338M gate-keeper mutation and maintain its decent kinase binding efficiency. However, a closer comparison of the kinase-binding activity of RL-45 with that of FDA-approved kinase inhibitor dasatinib shows that binding potency of RL-45 to the wild-type kinase inhibitor remains several order of magnitude weaker (IC50=0.021 ± 0.005 μM) compared to that of dasatinib (IC50=0.0004 ± 0.0002 μM). In the current work, we wanted to redesign potential kinase inhibitors, which can serve two purposes: a) similar to RL-45, these small molecules should be able to withstand drug resistance due to T338M mutation and b) these small molecules should have stronger kinase binding affinity than RL-45. As for assessing the ability of these small molecules to circumvent the drug resistance due to T338M mutation, Mondal et al14 had earlier shown that the central aromatic ring of RL-45 displays immaculate flexibility in bypassing any steric interaction resulting due to T338M mutation. Accordingly, in the current work, our fragmentbased drug design (FBDD) approach aims at keeping the central core of the RL-45 intact (see Figure 1), while introducing subtle changes in the various functional groups at different locations in the RL-45 core. We would demonstrate that the two derivatives of RL-45 (Figure 2), which we eventually have redesigned based on FBDD followed by QSAR and their varied chemistry (see Methods and figure S1 for details), are found to share a number of common protein-ligand interactions with RL-45, while bringing up additional interactions deemed important for enhancing kinase-ligand affinity than that in RL-45.

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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: The chemical structure of parent kinase inhibitor RL-45 and the two redesigned ligands LIG-1028 and LIG-1346. The dotted lines notify the position of introduction of new fragments around parent ligand RL-45

A close inspection of long equilibrium molecular dynamics simulation of RL-45 complexed with the src-kinase shows that the ligand remains localized in the cavity during the period of the simulation. Towards this end, we had individually parameterized the force fields of RL45 and the two newly redesigned derivatives parameters against ab-initio data25. Further Molecular Mechanics Poisson-Boltzman Surface Area (MM/PBSA)37 approach was employed to calculate the binding free energy in the present study. g_mmpbsa tool as developed by Kumari et al38 was used for this MM/PBSA calculation which combines Gromacs39 and APBS electrostatics program40 to calculate the total binding interaction energy between ligand and protein. The favorable binding free energies of each ligand in both wild type and mutated src kinase (SI Table 1), as obtained by MM/PBSA, suggested that the ligands are grossly stable inside the cavity. One of the main reasons for considering RL-45 ligand as our template is because of its proven ability7 in bypassing the resistance due to T338M gate-keeper residue mutation. Accordingly, we wanted to first confirm whether the two candidate RL-45 derivatives, which have been curated based on FBDD and subsequent QSAR-based screenings, maintain the 12 ACS Paragon Plus Environment

Page 12 of 25

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

The Journal of Physical Chemistry

ability to circumvent the T338M mutation. Accordingly, we individually computed the change in relative binding free energies (ΔΔG) of RL-45 and its two derivatives with wild type versus T338M mutated kinase using FEP/REST method (see Figure 3 for the thermodynamic cycle and computed ΔΔG values). This method has earlier found to provide superior free energy estimates by ensuring exhaustive sampling.14,19 Comparison of change in binding free energy of ligand molecules to the wild type versus mutated showed how RL-45 derivatives respond to T338M mutation. The table presented in Figure 3 enumerates the relative change of binding free-energy values (ΔΔG) obtained from FEP/REST simulations. The simulated ΔΔG value of -1 kcal/mol for RL-45 is in reasonable agreement with experimental value of 0.30 kcal/mol7,14. The 𝛥𝛥𝐺 values obtained from FEP/REST simulations showed us that all RL-45 derivatives could withstand the T338M mutation without much change in free energy.

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 3: Schematic of thermodynamic cycle used for the FEP/REST computations for the parent ligand RL-45 and the two newly redesigned ligands LIG-1028 and LIG-1346. Also shown in the table are the ΔΔG values for each of the three ligands, as obtained from FEP/REST calculation.

Mondal et al had earlier14 demonstrated that the flexibility of central phenyl ring of RL-45 is the key to maintain its similar binding ability to wild-type and T338M mutated kinase. Inspection of trajectories for all two RL-45 derivatives also reflected that the flexibility of phenyl ring has been retained in all the RL-45 derivative inhibitors. This is evident in the broad distributions of torsion angle of central phenyl ring of all RL-45 derivatives in both wild-type (𝜆 = 0) and mutated (𝜆 = 1) kinase (Figure 4). The observation of flexibility of central phenyl rings in all simulation related to RL-45 derivatives also lends credence to our efforts to individually parameterize the force fields small molecules against ab-initio data.25

Figure 4: Quantification of flexibility of phenyl ring of the ligand proximal to the gatekeeper residue studied using FEP/REST: The distribution of the torsional angle of the phenyl ring of ligand at 𝜆 =0.0 (wild-type) and

𝜆 =1.0 (T338M mutated) stage. The wide distribution of torsional angle in both the wild-type and mutant quantifies the high flexibility of the phenyl ring in all three ligands, which helps each of them to withstand the gatekeeper residue mutation.

The pair correlation between gatekeeper residue of kinase and RL-45 derivatives yielded similar results as in RL-45 (see Figure S2 in SI). Specifically, we observed the increase in the intensity of pair-correlation along with a shift in the peak towards smaller distance, which 14 ACS Paragon Plus Environment

Page 14 of 25

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

The Journal of Physical Chemistry

indicates lesser steric repulsion between methionine and ligands, as phenyl ring gets easily accommodated with methionine side-chain due to its flexibility. The major focus of the current work is to investigate whether these newly designed ligands can show relatively stronger binding affinity with the kinase compared to that of RL-45. However, considering the difficulty associated with computational estimation of absolute binding free energy of these ligands with the protein, we have taken a practical route where enhanced sampling MD simulations can be used to relatively rank these molecules on the basis of their resilience against a biasing potential. Towards this end, we follow the approach recently developed by Clark and coworkers17 where the metadynamics simulation had been employed to bias the ligand exit from the cavity and the ‘metadynamics simulation time’ needed for the ligand exit was used for ranking the ligand against their affinity to the pocket. Accordingly, in the current work, we have performed multiple independent metadynamics simulations exploring the unbinding pathway of our designed ligands through fragment-based drug discovery. We have taken the unbinding pathway of RL-45 as positive control for the interacting residues and intended to judge the designed ligands based on those interactions. Our results (see Figure 5) suggest that the fragments added influence the ligand recognition and expulsion almost in a similar way (with many common interacting residues as of RL-45), while the newly formed interactions tune the resilience of the ligand with the cavity.

15 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 5: Comparison of time profile of kinase-ligand distance for the parent ligand (RL-45) with the two redesigned ligand (LIG-1346 and LIG-1028), as obtained from metadynamics simulation for assessing the resilience of ligand against unbinding from its native binding pose (see method in SI for description of two collective variables used for metadynamics simulation). For each ligand, three independent metadynamics simulation trajectories are shown. The ligand unbinding time (characterized by considering 2 nm cutoff distance) for each ligand is shown in the table included in the figure. The longer time for ligand exit in case of LIG-1028 and LIG-1346 relative to RL-45 suggests newly redesigned ligands will be more resilient against unbinding from kinase.

The table presented inside Figure 5 compares the “metadynamics simulation time” needed for each of the three ligands to exit the kinase cavity, as obtained multiple independent trajectories. Previously, the metadynamics simulation time required to unbind the ligand from its native binding pose has been treated as a decent proxy of the true ligand-unbinding time (which otherwise needs very slow metadynamics biasing33,34) and have been used as an efficient metric to screen the small molecules by Clark et al.17 By comparison of simulated ‘metadynamics time’, a significant increase in unbinding times for LIG-1028 and LIG-1346

16 ACS Paragon Plus Environment

Page 16 of 25

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

The Journal of Physical Chemistry

from the binding pocket of the kinase relative to RL-45 is observed, which signifies an increased binding affinity of these two ligands towards cSrc kinase over RL-45. Together, the screening based upon metadynamics simulation predicts LIG-1028 and LIG-1346 to have potentially superior kinase-binding ability than RL-45. We emphasize that the ‘metadynamics run time’ reported here are not true ligand-residence time. Similar to Clark et al17, we have also suggested in this work that this ‘metadynamics run time’ can be practical yardstick for screening the ligands. However we are also aware of the fact that the simulations with similar runtime could have different placement of bias potential. Hence, to check whether this comparison via ‘metadynamics runtime’ translates to stronger binding affinity, we have computed approximate free energy profiles of binding affinity of these ligands from these metadynamics simulation trajectories by adding up the deposited bias. The comparison of free energy profiles, as shown in Figure S3 in supporting information, show trends consistent with the respective metadynamics run time, the one used as the practical metric for screening the ligand.

To assess whether the predicted trend of ligand potency remains robust irrespective of bias deposition interval, we have additionally performed metadynamics simulation for two ligands (RL-45 and LIG-1346) at 2.5 times longer bias deposition interval than the original one (i.e. 2 ps contrary to 0.8 ps). The ‘metadynamics run time’ among two ligands at slower bias frequency is compared in figure S4 in supporting information. As expected, while the ‘metadynamics run time’ for unbinding the ligand increases with longer bias interval for both the ligands, the originally predicted trend of superior binding potency of the designed ligand LIG-1346 over RL-45 holds true in both bias deposition interval. The metadynamics simulation trajectories render key molecular insights on the higher resilience of the two screened potential inhibitors, namely LIG-1346 and LIG-1028 against

17 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

unbinding, compared to the parent inhibitor RL-45. As a set of common interactions in all three ligands, the central phenyl ring, as observed from the bound pose of all candidate inhibitors, interacted with Leu393, Val325, alkyl chain of Lys295, Thr338 and Val323 of active site, making a stable hydrophobic sub-pocket. However, in addition to these favorable interactions common with protein-RL-45 complex, the two newly redesigned ligand LIG1346 and LIG-1028 reinforce more consistent hydrogen bonds with the pocket residues and establish a key hydrophobic interaction during the period of metadynamics driven expulsion process. Specifically, as depicted in figure 6, the two newly redesigned ligands maintain hydrogen bonds with Met341, Glu310 and Asp404 for significantly longer times than the parent inhibitor RL-45. Most importantly, the keto and primary amine group of newly added 4,4-dimethylisoquinoline-1,3(2H,4H)-dione fragment of LIG-1346 makes additional resilient hydrogen bonds with Gly274 (see Figure 6), which plays an important role in holding the ligand against unbinding from the pocket. We have also observed occasional hydrogen bonds between Gln275 of the protein with LIG-1346. Additionally the 4,4-dimethylisoquinoline1,3(2H,4H)-dione fragment makes 𝜋 ― 𝜋 stacking with Phe405, which contributes to LIG1346’s stronger affinity to the active site. Similarly, just before the ligand expulsion from active site, the newly introduced fragment in LIG-1028 made intermittent hbonds with Glu310, Lys295, Ala306, otherwise absent in RL-45-protein complex. More significantly, this newly introduced fragment, during its unbinding pathway, showed binding pocket resilience through a shifted hydrophobic paradigm where it was observed to form hydrophobic contacts through Leu322, Val323, Met314, Phe405, and Ala403 with pyrrolidine fragment. Together, these stabilizing interactions present in these redesigned small-molecules are projected to enhance the protein-ligand affinity.

18 ACS Paragon Plus Environment

Page 18 of 25

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

The Journal of Physical Chemistry

Figure 6: The snapshots capturing the key ligand interactions for each of the three ligands. The time profile of crucial hydrogen bonds for all three ligands is displayed. Also shown is the hydrophobic pocket in proteinLIG1028 complex. The interacting functional groups of LIG-1028, which participate in the shifted hydrophobic paradigm is shown in magenta.

Conclusion: It is widely recognized that drug discovery pertaining to kinases, in general, is a time consuming, expensive and extremely complex process. There are ample interests in probing unknown ‘chemical space’ to identify novel bioactive compounds with diverse scaffolds to improve the effectiveness of drug discovery. However, the field of computational drug design often faces challenges with accurately and reliably predicting the resilience of protein-ligand bound complex. The uncertainty associated with this final step significantly reduces the utility of docking calculations in actual projects, where actionable drug discovery efforts such as synthetic prioritization require a high confidence in the predicted binding 19 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

mode of the ligand. However, with the advancements in computing hardwares, recent years have seen the emergence of unprecedently long molecular dynamics simulations and enhanced sampling techniques, which have been able to successfully capture the full process of ligand entry34,41–43 to and ligand exit34,44–46 from the protein-binding cavity, including that in kinases.34,41 As a possible middle ground, recent computational drug discovery protocols are increasingly introducing explicit solvent molecular dynamics simulations to discriminate between the relatively small number of top ranked poses produced by the docking projects.17,47 The current work cross-validates the resilience of the src kinase-ligand complex developed by fragment-based drug design and in doing so, takes a step forward by integrating the explicit solvent MD simulation with a pair of free-energy based enhanced sampling techniques, namely metadynamics and FEP/REST. While this work is inspired form our previous research initiative14, the sole purpose of our early work was to provide an mechanistic interpretation of the two popular ligands’ (dasatinib and RL-45) different capability at circumventing T338M mutation. The current work takes an important step forward and attempts to propose, via combination of fragment based drug design and enhanced sampling techniques, new kinase inhibitor derivatives which can retain the ability of one of ligand’s gate-keeper residue mutation circumvention ability plus improve the kinase binding potency. Specifically, we show how to redesign a previously developed kinase inhibitor RL-45 by optimizing both its ability of withstanding drug-resistance related to gatekeeper mutation and its resilience against stronger binding affinity towards kinase. The newly designed kinase inhibitors retain the flexibility of parent drug, which is crucial for circumventing gate-keeper residue mutation. Moreover, these two newly designed compounds, via long-standing hydrogen bonds, specifically involving Met341 and Asp404, can assert strong affinity with the Src-kinase. In addition to that, exploitation of the residue

20 ACS Paragon Plus Environment

Page 20 of 25

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

The Journal of Physical Chemistry

hydrophobicity inside the binding pocket, as noticed in the case of LIG-1028, can prove fruitful in retaining potential inhibitors inside binding pocket. Finally we note that the ideal validation of the current computational results would have been experimental demonstration of results presented in this work. However, due to academic limitations of the authors, it is not possible at present to conduct any wet lab experiments to validate the hypothesis. However, we hope that our results, once published, would motivate medicinal chemists to synthesize the designed ligands and assess the potency of these new ligands.

Supplemental Information: Chemical structures of the sort-listed molecules from QSAR, Table with MM/PBSA results and figures with pair correlation function of gate keeper residue with each of the ligand, free energy profiles and ligand unbinding time profile at slower bias frequency (PDF) QSAR dataset in sdf format (zip file) Acknowledgement: This work was supported by computing resources obtained from shared facility of TIFR Hyderabad, India. JM would like to acknowledge research intramural research grants obtained from TIFR, India, Ramanujan Fellowship and Early Career Research funds provided by the Department of Science and Technology (DST) of India (ECR/2016/000672). References (1)

Treiber, D. K.; Shah, N. P. Ins and Outs of Kinase DFG Motifs. Chem. Biol. 2013, 20 (6), 745– 746. https://doi.org/10.1016/j.chembiol.2013.06.001.

(2)

Hari, S.; Merritt, E.; Maly, D. Sequence Determinants of a Specific Inactive Protein Kinase Conformation. Chem. Biol. 2013, 20 (6), 806–815. https://doi.org/10.1016/j.chembiol.2013.05.005.

(3)

Cohen, P. Protein Kinases--the Major Drug Targets of the Twenty-First Century? Nat. Rev.

21 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

Drug Discov. 2002, 1 (4), 309–315. https://doi.org/10.1038/nrd773. (4)

Zhang, J.; Yang, P. L.; Gray, N. S. Targeting Cancer with Small Molecule Kinase Inhibitors. Nat. Rev. Cancer 2009, 9 (1), 28–39. https://doi.org/10.1038/nrc2559.

(5)

Backes, A. C.; Zech, B.; Felber, B.; Klebl, B.; Muller, G. Small-Molecule Inhibitors Binding to Protein Kinases. Part I: Exceptions from the Traditional Pharmacophore Approach of Type I Inhibition. Expert Opin. Drug Discov. 2008, 3 (12), 1409–1425. https://doi.org/10.1517/17460440802579975.

(6)

Bikker, J. A.; Brooijmans, N.; Wissner, A.; Mansour, T. S. Kinase Domain Mutations in Cancer : Implications for Small Molecule Drug Design Strategies. J.Med.Chem. 2009, 52 (6), 1493– 1509.

(7)

Getlik, M.; Grütter, C.; Simard, J. R.; Klüter, S.; Rabiller, M.; Rode, H. B.; Robubi, A.; Rauh, D. Hybrid Compound Design to Overcome the Gatekeeper T338M Mutation in CSrc. J. Med. Chem. 2009, 52 (13), 3915–3926. https://doi.org/10.1021/jm9002928.

(8)

Daub, H.; Specht, K.; Ullrich, A. Strategies to Overcome Resistance to Targeted Protein Kinase Inhibitors. Nat. Rev. Drug Discov. 2004, 3 (12), 1001–1010. https://doi.org/10.1038/nrd1579.

(9)

Carter, T. A.; Wodicka, L. M.; Shah, N. P.; Velasco, A. M.; Fabian, M. A.; Treiber, D. K.; Milanov, Z. V; Atteridge, C. E.; Iii, W. H. B.; Edeen, P. T.; et al. Inhibition of Drug-Resistant Mutants of ABL , KIT , and EGF Receptor Kinases. Proc Natl Acad Sci U S A 2005, 102 (31), 11011–11016.

(10)

Gibbons, D. L.; Pricl, S.; Kantarjian, H.; Cortes, J.; Quintas-Cardama, A. The Rise and Fall of Gatekeeper Mutations? The BCR-ABL1 T315I Paradigm. Cancer 2012, 118 (2), 293–299. https://doi.org/10.1002/cncr.26225.

(11)

Young, M. A.; Shah, N. P.; Chao, L. H.; Seeliger, M.; Milanov, Z. V; Biggs, W. H.; Treiber, D. K.; Patel, H. K.; Zarrinkar, P. P.; Lockhart, D. J.; et al. Structure of the Kinase Domain of an Imatinib-Resistant Abl Mutant in Complex with the Aurora Kinase Inhibitor VX-680. Cancer Res. 2006, 66 (2), 1007–1014. https://doi.org/10.1158/0008-5472.CAN-05-2788.

(12)

Karaman, M. W.; Herrgard, S.; Treiber, D. K.; Gallant, P.; Atteridge, C. E.; Campbell, B. T.; Chan, K. W.; Ciceri, P.; Davis, M. I.; Edeen, P. T.; et al. A Quantitative Analysis of Kinase Inhibitor Selectivity. Nat. Biotechnol. 2008, 26, 127.

(13)

Mayeda, M.; Weisberg, E.; Moreno, D.; Khosravi-Far, R.; Melo, J. V.; Debiec-Rychter, M.; Ray, A.; Gray, N.; Azam, M.; Manley, P. W.; et al. Discovery of a Small-Molecule Type II Inhibitor of Wild-Type and Gatekeeper Mutants of BCR-ABL, PDGFR , Kit, and Src Kinases: Novel Type II Inhibitor of Gatekeeper Mutants. Blood 2010, 115 (21), 4206–4216. https://doi.org/10.1182/blood-2009-11-251751.

(14)

Mondal, J.; Tiwary, P.; Berne, B. J. How a Kinase Inhibitor Withstands Gatekeeper Residue Mutations. J. Am. Chem. Soc. 2016, 138 (13), 4608–4615. https://doi.org/10.1021/jacs.6b01232.

(15)

Doak, B. C.; Norton, R. S.; Scanlon, M. J. The Ways and Means of Fragment-Based Drug Design. Pharmacol. Ther. 2016, 167, 28–37. https://doi.org/https://doi.org/10.1016/j.pharmthera.2016.07.003.

(16)

Ajmani, S.; Jadhav, K.; Kulkarni, S. A. Three-Dimensional QSAR Using the k-Nearest Neighbor Method and Its Interpretation. J. Chem. Inf. Model. 2006, 46 (1), 24–31. https://doi.org/10.1021/ci0501286.

22 ACS Paragon Plus Environment

Page 22 of 25

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

The Journal of Physical Chemistry

(17)

Clark, A. J.; Tiwary, P.; Borrelli, K.; Feng, S.; Miller, E. B.; Abel, R.; Friesner, R. A.; Berne, B. J. Prediction of Protein–Ligand Binding Poses via a Combination of Induced Fit Docking and Metadynamics Simulations. J. Chem. Theory Comput. 2016, 12 (6), 2990–2998. https://doi.org/10.1021/acs.jctc.6b00201.

(18)

Wang, L.; Berne, B. J.; Friesner, R. A. On Achieving High Accuracy and Reliability in the Calculation of Relative Protein-Ligand Binding Affinities. Proc. Natl. Acad. Sci. 2012, 109 (6), 1937–1942.

(19)

Wang, L.; Wu, Y.; Deng, Y.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.; Dahlgren, M. K.; Greenwood, J.; et al. Accurate and Reliable Prediction of Relative Ligand Binding Potency in Prospective Drug Discovery by Way of a Modern Free-Energy Calculation Protocol and Force Field. J. Am. Chem. Soc. 2015, 137, 2695–2703.

(20)

Yang, G.-F.; Jiang, W.; Zhu, X.-L.; Hao, G.-F.; Guo, F.-B.; Wu, F.-X.; Ye, Y.-N. ACFIS: A Web Server for Fragment-Based Drug Discovery. Nucleic Acids Res. 2016, 44 (W1), W550–W556. https://doi.org/10.1093/nar/gkw393.

(21)

Gaulton, A.; Bellis, L. J.; Bento, A. P.; Chambers, J.; Davies, M.; Hersey, A.; Light, Y.; McGlinchey, S.; Michalovich, D.; Al-Lazikani, B.; et al. ChEMBL: A Large-Scale Bioactivity Database for Drug Discovery. Nucleic Acids Res. 2012, 40 (Database issue), D1100–D1107. https://doi.org/10.1093/nar/gkr777.

(22)

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

(23)

Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E. M.; Mittal, J.; Feig, M.; MacKerell, A. D. Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain Χ1 and Χ2 Dihedral Angles. J.Chem.Theory Comput. 2012, 8 (9), 3257–3273. https://doi.org/10.1021/ct300400x.

(24)

Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM General Force Field: A Force Field for Drug-like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31 (4), 671–690. https://doi.org/10.1002/jcc.21367.

(25)

Huang, L.; Roux, B. Automated Force Field Parameterization for Nonpolarizable and Polarizable Atomic Models Based on Ab Initio Target Data. J. Chem. Theory Comput. 2013, 9 (8), 3543–3556. https://doi.org/10.1021/ct4003477.

(26)

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. W. Comparison of Simple Potential Functions for Simulating Liquid Water. J.Chem.Phys. 1983, 79, 926–935.

(27)

Limongelli, V.; Bonomi, M.; Marinelli, L.; Gervasio, F. L.; Cavalli, A.; Novellino, E.; Parrinello, M. Molecular Basis of Cyclooxygenase Enzymes (COXs) Selective Inhibition. Proc. Natl. Acad. Sci. U. S. A. 2010, 107 (12), 5411–5416. https://doi.org/10.1073/pnas.0913377107.

(28)

Gervasio, F. L.; Laio, A.; Parrinello, M. Flexible Docking in Solution Using Metadynamics. J. Am. Chem. Soc. 2005, 127 (8), 2600–2607. https://doi.org/10.1021/ja0445950.

(29)

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

23 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

(30)

Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; et al. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comp. Phys. Comm. 2009, 180, 1961–1972.

(31)

Copeland, R. A.; Pompliano, D. L.; Meek, T. A. Drug-Target Residence Time and Its Implications for Lead Optimization. Nat. Rev. Drug Discov. 2006, 5, 730–739.

(32)

Tiwary, P.; Parrinello, M. From Metadynamics to Dynamics. Phys Rev Lett 2013, 111 (23), 230602–230606. https://doi.org/10.1103/PhysRevLett.111.230602.

(33)

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

(34)

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

(35)

Salvalaglio, M.; Tiwary, P.; Parrinello, M. Assessing the Reliability of the Dynamics Reconstructed from Metadynamics. J.Chem.Theory Comput. 2014, 10 (4), 1420–1425.

(36)

Gapsys, V.; Michielssens, S.; Seeliger, D.; de Groot, B. L. Pmx: Automated Protein Structure and Topology Generation for Alchemical Perturbations. J. Comput. Chem. 2015, 36 (5), 348– 354. https://doi.org/10.1002/jcc.23804.

(37)

Genheden, S.; Ryde, U. The MM/PBSA and MM/GBSA Methods to Estimate Ligand-Binding Affinities. Expert Opin. Drug Discov. 2015, 10 (5), 449–461.

(38)

Kumari, R.; Kumar, R.; Lynn, A. G_mmpbsa—A GROMACS Tool for High-Throughput MMPBSA Calculations. J. Chem. Inf. Model. 2014, 54 (7), 1951–1962. https://doi.org/10.1021/ci500020m.

(39)

Hess, B.; Kutzner, C.; van der Spoel, D. GROMACS 4: Algorithms for Highly Efficient, LoadBalanced, and Scalable Molecular Simulation. J. Chem. Theor. Comput. 2008, 4, 435–447.

(40)

Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Electrostatics of Nanosystems: Application to Microtubules and the Ribosome. Proc. Natl. Acad. Sci. 2001, 98 (18), 10037 LP10041. https://doi.org/10.1073/pnas.181342398.

(41)

Shan, Y.; Kim, E. T.; Eastwood, M. P.; Dror, R. O.; Seeliger, M. A.; Shaw, D. E. How Does a Drug Molecule Find Its Target Binding Site? J. Am. Chem. Soc. 2011, 133, 9181–9183.

(42)

Mondal, J.; Ahalawat, N.; Pandit, S.; Kay, L. E.; Vallurupalli, P. Atomic Resolution Mechanism of Ligand Binding to a Solvent Inaccessible Cavity in T4 Lysozyme. PLOS Comput. Biol. 2018, 14 (5), e1006180.

(43)

Ahalawat, N.; Mondal, J. Mapping the Substrate Recognition Pathway in Cytochrome P450. J. Am. Chem. Soc. 2018, 140 (50), 17743–17752. https://doi.org/10.1021/jacs.8b10840.

(44)

Feher, V. A.; Schiffer, J. M.; Mermelstein, D. J.; Mih, N.; Pierce, L. C. T.; McCammon, J. A.; Amaro, R. E. Mechanisms for Benzene Dissociation through the Excited State of T4 Lysozyme L99A Mutant. Biophys. J. 2019, 116 (2), 205–214. https://doi.org/10.1016/j.bpj.2018.09.035.

(45)

Lotz, S. D.; Dickson, A. Unbiased Molecular Dynamics of 11 Min Timescale Drug Unbinding Reveals Transition State Stabilizing Interactions. J. Am. Chem. Soc. 2018, 140, 618–628. https://doi.org/10.1021/jacs.7b08572.

24 ACS Paragon Plus Environment

Page 24 of 25

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

The Journal of Physical Chemistry

(46)

Nunes-Alves, A.; Zuckerman, D. M.; Arantes, G. M. Escape of a Small Molecule from Inside T4 Lysozyme by Multiple Pathways. Biophys. J. 2018, 114 (5), 1058–1066. https://doi.org/10.1016/j.bpj.2018.01.014.

(47)

Śledź, P.; Caflisch, A. Protein Structure-Based Drug Design: From Docking to Molecular Dynamics. Curr. Opin. Struct. Biol. 2018, 48, 93–102. https://doi.org/10.1016/j.sbi.2017.10.010.

TOC Graphics

25 ACS Paragon Plus Environment