Towards Atomistic Modelling of Irreversible Covalent Inhibitor Binding

Subscriber access provided by Bibliothèque de l'Université Paris-Sud

Pharmaceutical Modeling

Towards Atomistic Modelling of Irreversible Covalent Inhibitor Binding Kinetics Haoyu S. Yu, Cen Gao, Dmitry Lupyan, Yujie Wu, Takayuki Kimura, Chuanjie Wu, Leif jacobson, Edward D. Harder, Robert Abel, and Lingle Wang J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.9b00268 • Publication Date (Web): 19 Aug 2019 Downloaded from pubs.acs.org on August 20, 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 39 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

Towards Atomistic Modelling of Irreversible Covalent Inhibitor Binding Kinetics Haoyu S. Yu,a, Cen Gao, b Dmitry Lupyan,a Yujie Wu,a Takayuki Kimura,c Chuanjie Wu, a Leif Jacobson,a Edward Harder, a Robert Abel,a and Lingle Wang*a

a) Schrodinger, Inc., 120 West 45th Street, New York, New York 10036, United States b) Eli Lilly and Company, Lilly Corporate Center, Indianapolis, IN 46285, United States c) Schrodinger, Inc., 101 SW Main Street, Suite 1300, Portland, Oregon, 97204, United States Corresponding Author *E-mail:[email protected] Abstract Covalent Inhibitors have emerged as an important drug class in recent years, largely due to their many unique advantages as compared to noncovalent inhibitors, including longer duration of action, lower prolonged systemic exposure, higher potency and selectivity. However, the potential off-target toxicity of covalent inhibitors, particularly of irreversible covalent inhibitors, represents a great challenge in covalent drug development. Therefore, accurate calculation of protein covalent inhibitor reaction kinetics to guide the design of selective inhibitors would greatly benefit covalent drug discovery efforts. In the present paper, we present a computational method to calculate the relative reaction kinetics between congeneric irreversible covalent inhibitors and their protein receptors. The method combines density functional theory calculations of the transition state barrier height of the rate limiting step for reaction between the warhead of the inhibitor and a single protein residue, and molecular mechanical-based free energy calculations to account for the interactions between the ligand in transition state and the protein environment. The method was tested on 4 pharmaceutically interesting irreversible covalent binding systems

ACS Paragon Plus Environment

1

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

Page 2 of 39

involving 28 ligands; the mean unsigned error (MUE) of the relative reaction rate for all pairs of ligands between the predictions and experimental results for these tested systems is 0.79 log units. This is to our knowledge the first time where the reaction kinetics of protein irreversible covalent inhibition have been directly calculated with physics-based free energy calculation methods and transition state theory. We anticipate the outstanding accuracy demonstrated here across a broad range of target classes will have a strong impact on the design of selective covalent inhibitors.


2

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


1. Introduction Irreversible covalent inhibitors represent a distinct and important class of pharmaceutical drugs. A large fraction of Food and Drug Administration (FDA) approved drugs achieve their therapeutic effects through formation of a covalent bond with a target protein receptor. Notable examples of available treatments utilizing this mode of action include the aspirin and penicillin,1 the recent ‘blockbuster’ proton pump inhibitor omeprazole and antiplatelet drug clopidogrel,2 and the most recently developed kinase inhibitors afatinib,3 ibrutinib4 and rociletinib5 for cancer therapies. Despite the large number of successful covalent drugs, drug designers remain reluctant to pursue a covalent mode of action in drug discovery programs, largely owing to perceived increased risk of off-target toxicity. This concern is not unfounded; pioneering work conducted in the 1970s regarding the hepatotoxic properties of drug-like fragments including bromobenzene and acetaminophen provided solid evidence of off-target toxicity for compounds with chemically reactive metabolites.6 Due to these concerns, essentially all earlier covalent drugs were discovered through screening in biological assays, with their covalent mechanism of action established afterwards. In recent years, covalent inhibitors have experienced a resurgence of interest, mainly due to their unique advantages as compared to non-covalent ligands. These unique advantages include longer duration of action, lower prolonged systemic exposure, and higher potency and selectivity.7 A notable recent success story has been the design of targeted covalent kinase inhibitors (TCKI). The warhead of TCKI is carefully tuned to achieve low reactivity with other off-target proteins while the remainder of the ligand forms highly complementary non-covalent interactions with the target protein binding pocket. In this way, TCKI can achieve favourable selectivity profile to avoid off-target toxicity. Notable examples of TCKIs include epidermal growth factor receptor (EGFR) inhibitors afatinib3 and rociletinib,5 and Bruton tyrosine kinase (BTK) inhibitor ibrutinib4 for cancer therapies. Depending on the reaction barriers covalent inhibitors can bind either reversibly or irreversibly. For reversible covalent inhibitors, the reaction barrier for the formation and


3


Page 4 of 39

breaking of the chemical bond between the ligand warhead and targeted protein residue is sufficiently low such that the formation of the covalent linkage is reversible on the time scale of an experimental assay. For such inhibitors, the covalent complex and the protein and ligand separated in solution are in equilibrium and the binding free energy is a useful quantity to evaluate the efficacy of the ligand. On the other hand, for irreversible covalent inhibitors, the reaction barrier for the breaking of the covalent bond between the ligand warhead and the targeted protein residue is sufficiently high to make the covalent bond formation reaction effectively irreversible on the time scale of an experimental assay. For such inhibitors, their binding affinities are time-dependent, and the efficacies are determined by kinetics. The quantities used to describe the kinetics of the binding process are the rate of the covalent bond formation ( kinact ) as well as the forward ( kon ) and reverse rate ( koff ) of non-covalent binding of the inhibitor prior to formation of the covalent linkage. The desired targeted covalent inhibitor (TCI) should have slow kinact to avoid off-target toxicity while manifesting strong binding affinity for the pre-reaction non-covalent complex. In a previous study, we have designed a method to calculate the relative binding free energies between congeneric reversible covalent inhibitors,8 and the same method was also used by Chatterjee et al to study the selectivity of reversible covalent inhibitors binding to two highly homologous proteases, calpain-1 and calpain-2.9 In the present paper, we extend the method to the calculation of the relative binding kinetics of irreversible covalent inhibitors. The general strategy is to use DFT to model the transition state of a model reaction between the ligand warhead and a single protein residue it forms a covalent bond with, and then to use the existing free energy perturbation method to model the interactions between the protein environment and the ligands. The treatment avoids the need of the hybrid quantum mechanics/molecular mechanics approach (QM/MM) to model the reaction between the ligand and the entire protein. In particular, the method first calculates the transition state barrier height of the rate limiting step for reaction between the warhead of the inhibitor and the targeted protein residue with which the warhead will react, and then calculates the effect of protein environment on the


4

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


transition state barrier through molecular mechanics-based free energy calculations. The method was tested on 4 pharmaceutically interesting irreversible covalent binding systems involving 28 ligands and the calculated relative binding kinetics agree well with experimental data, with the MUE of Δlog(kinact/Ki) (the relative reaction rate between a pair of congeneric ligands in log unit) between the predicted and experimental results of 0.75 log unit. This is to our knowledge the first time where the reaction kinetics of protein irreversible covalent inhibition has been directly calculated with physics-based methods. 2. Methods 2.1 Theory The binding process of an irreversible covalent ligand to the protein target involves two steps, formation of a non-covalently bound complex and formation of a covalent bond. Covalent inhibitor design requires careful optimization of the non-covalent binding affinity

Ki and the reactivity of the warhead kinact . In Figure 1, we show the potential energy curve for irreversible covalent ligand binding to the target protein and the corresponding binding kinetics. Starting from the unbound state (P+L), van der Waals and electrostatic interactions drive the formation of the protein-ligand non-covalent complex (P.L) with binding free energy of Gi  RT ln Ki , where Ki is the equilibrium constant, which is equal to the ratio of the unbinding and binding rate constant, koff / kon . The non-covalent complex (P.L) will further go through a transition state (P-L≠) with reaction barrier of G  to form the protein-ligand covalent complex (P-L). According to the transition state theory, the rate of covalent bond formation, which is also the rate the enzyme (protein) becomes inactive, kinact , is dependent on the reaction barrier through 

 [P  L] kinact  Ae G / RT  A [P.L]

(1)

where A is a constant value from transition state theory, ΔG≠ is standard Gibbs free energy of activation10, and [P-L] and [P.L] stand for the equilibration concentrations of the protein-ligand complex in the transition state and in the non-covalent bound form. The apparent, or overall, rate for the formation of the protein-ligand covalent complex from


5


Page 6 of 39

separated protein and ligand in solution, kinact / Ki , is the key property to optimize in irreversible covalent ligand design. The ideal ligand should have slow kinact , to avoid offtarget related toxicities, and favorable non-covalent binding affinity to optimize the selectivity. From Figure 1, it is clear that the overall reaction rate kinact / Ki is related to the difference between the free energy of activation of the transition state and free energy of the reactant, ie,

kinact / Ki  Ae

(G  Gi )/RT

 Ae

 / RT Gapp

A

[P  L] [P][L]

(2)

  G   Gi stands for the apparent free energy of activation which is the where Gapp

free energy difference between the transition state10 and the reactant. In equation 2, [P-L], [P], and [L1] stand for equilibrium concentrations of protein-ligand complex in transition state, protein, and ligand respectively. Accurate calculation of irreversible covalent binding kinetics requires the modeling of the transition state of the protein-ligand covalent bond formation reaction, which is very challenging. In the following, we propose an approximate method to calculate the relative binding kinetics between two structurally similar irreversible covalent ligands. The proposed method combines density functional theory (DFT) calculations of the transition state barrier height for the reaction between the inhibitor and a “proxy” single protein residue (or a fragment of the single residue) with which it forms covalent bond, and molecular mechanics-based free energy calculations to account for the effect of protein environment on the transition state barrier. There are a few assumptions in this model: (1) the rate limiting step in solution is the same as the one in protein; (2) the geometry of the transition state is not significantly altered by the presence of the protein; (3) only the electrostatic and van der Waals interactions between the protein and the ligand in the transition state geometry change the reaction barrier height in the protein binding pocket; (4) when the perturbation region is far from the warhead region of the covalent ligand, for example, separated for more than 6 carbons and the perturbed region is not conjugated


6

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


with the warhead, we further assume that the transition state geometry will not vary with respect to different R groups on the tail of the ligand. Consider the relative binding kinetics between two structurally related ligands. The ratio of the apparent rate of formation between the two ligands can be written as

[P  L1] [P][L1]

[R  L1] [R][L1]

[P  L1]

(2) (1) kinact [P  L2 ] kinact :    (3) Ki(1) Ki(2) [P  L2 ] [R  L2 ] [R  L1] [P][L2 ] [R][L2 ] [R  L ] 2

  The second step of the equation is obtained by inserting [R  L1] / ([R][R  L2 ] ) in both

the denominator and numerator, where R stands for any model of the single protein residue the ligand warhead reacts with to form a covalent bond, and [R-L1] and [R] stand for the equilibrium concentrations of modeled residue-ligand 1 complex in the transition state and the single protein residue itself, respectively. In Equation 3, the overall reaction rate is decomposed into the product of two terms. The first term is the ratio of equilibrium concentrations of the transition state and the reactants for the reaction between ligand and a small model of the single reactive protein residue amenable to DFT computations while the second term accounts for the difference between the entire protein and this model. In Figure 2, the free energy difference between the transition state of protein-ligand covalent complex and the reactant is decomposed into two parts, first part is the free energy difference between the transition state and its reactants for the single protein residue and the ligand reaction, the second part is a correction term that accounts for the different interactions with ligand in transition state between the entire protein and a single protein residue. In our protocol the two terms appearing in Equation 3 are computed by two separate calculations. In particular, the relative rates for the reactions between the two ligands with a single protein residue (the first term of Equation 3) is computed by DFT transition state calculations,11 while the second term in Equation 3 can be computed by free energy


7


perturbation (FEP)

methods.12

The numerator in the second term (

Page 8 of 39

[P  L1]

[P  L2 ]

) is related to

the free energy of transforming protein-ligand 1 covalent complex into protein-ligand 2 covalent complex in the transition state, where the denominator in the second term of Equation 3 (

[R  L1]

[R  L2 ]

) is related to the free energy to transform a single protein residue-

ligand 1 covalent complex to a single protein residue-ligand 2 covalent complex in the transition state. Therefore, by combining DFT transition state energy calculations and FEP calculations in the transition state, we can compute the relative binding kinetics between two structurally related ligands. To demonstrate how the proposed method works for relative binding kinetics calculations, a thermodynamic cycle connecting the transition states and reactants for a pair of irreversible covalent inhibitors (Ligand 36 and Ligand 42) binding to BTK protein is shown in Figure 3.13 In this case, the acrylamide of the ligand reacts with the cysteine residue of the protein through Michael addition, one of the most common reactions in covalent inhibitor designs.14 The rate limiting step of this type of reaction is believed to be reaction of acrylamide with the deprotonated form of cysteine.15 The fluorine atom on the fused aromatic moiety in ligand 36 is transformed to a chlorine atom in Ligand 42 and an additional fluorine atom is also introduced in the terminal phenyl ring in Ligand 42. In Figure 3, we decompose the overall energetics for the reaction between ligand and the protein into two components. The first component is the reaction between the ligands and a single cysteine residue in the solvent with the reaction barriers of ΔG1(36) and ΔG1(42) respectively, which can be computed by DFT. In the second step, we compute the relative contribution from the protein environment to the transition state energy barriers of the two ligands, denoted by ΔG2(36) and ΔG2(42) respectively in the figure. In particular, we have computed this quantity by transforming protein-ligand 36 complex into protein-ligand 42 complex while restraining the ligand-protein interaction geometry to the transition state, and likewise repeating the analogous calculation for protein residue 442/Ligand 36 complex as well as protein residue 442/Ligand 42 complex in transition states. Since free


8

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


energy is a state function, the overall relative free energy of activation between the two ligands can be computed by these two steps.   Gapp (36)  Gapp (42)  (G1(36)  G 2 (36))  (G1(42)  G 2 (42))

 (G1(36)  G1(42))  (Gcomplex  Gsolvent )

(4)

  G   Gi is defined in Equation 2. The The apparent free energy of activation Gapp

first parenthesis in the last step of Equation 4 can be computed by calculating transition state barriers with DFT, for example by using the AutoTS11 program, and the second parenthesis can be computed by free energy calculation, for example by using FEP+.16-18 It is worth mentioning that the transition state is located on the saddle point of the potential energy surface (PES), which is not a stable structure for molecular dynamics simulations. To model the protein-ligand complex and the single residue-ligand complex in the transition states, we modified the force field parameters and added additional restraints in the molecular mechanics potential energies to restrict the complex structure in a geometry corresponding to the transition states. This is consistent with the definition of Gibbs activation free energy in the transition state where the degree of freedom along the reaction coordination is frozen.10 While a more sophisticated approach utilizing special force field to model the transition state has been attempted in the past,19 we found that modified equilibrium bond distances and angles in the force fields together with simple harmonic restraints on the key bond distances and angles in molecular dynamics (MD) simulations are sufficient to maintain the transition state geometry for the type of reactions modelled in this paper. In the case of the BTK system mentioned above, a harmonic restraint with force constant of 300 kcal/mol/Å2 was used to maintain the bond distance between the carbon atom on the ligand warhead and the sulphur atom from cysteine residue 442 that form a covalent bond in transition state in the FEP simulations. In the course of FEP simulations, we observed that this covalent bond distance remains essentially unchanged (fluctuation is less than 0.1 Å around the equilibrium distance of 1.9 Å) during the molecular dynamic trajectories. For the other types of reactions studied in this paper, including carbamylation reaction, the equilibrium bond distances, angles and torsions spanning the covalent bond between the protein and ligand were modified in the force


9


Page 10 of 39

fields to match the transition state geometry. In particular, the equilibrium distances of amide C-N bond in the ligand, the C-O bond between ligand and serine residue, and the NH bond between the ligand and the serine residue, were modified to 1.6 Å, 1.9 Å, and 1.2 Å respectively as shown in Figure 4. In addition, the equilibrium angles and torsions spanning the planar four-membered ring (C-N bond from the ligand and O-H bond from serine residue) were all modified in the force fields to match the transition state geometry. We observed that the protein-ligand complex remained in the transition state during the course of FEP simulations. It should be noted that the above model for calculating the irreversible covalent inhibitor binding kinetic has two assumptions. First, the reaction mechanism can be fully described by a single residue with which the ligand reacts to form a covalent bond. For the BTK systems in the section 3.1 and KRAS system in section 3.3, we performed QM/MM calculations to search for the transition state for the reaction between the ligand and the protein binding pocket, and the geometry of the transition state in QM/MM matches that from DFT calculations of the model reaction. Therefore, this assumption holds for the type of reactions we modelled in this paper. More information about this comparison can be found in the supporting information. In general, if the reaction between the ligand and the entire protein goes through a different transition state (a different mechanism) than that between the same ligand and a single protein residue, we would not expect our method to be able to predict the relative kinetics among irreversible covalent inhibitors. Second, the observed binding kinetics is assumed to be determined by a single rate limiting step of the entire chain of elementary reactions. If multiple steps of reactions contribute significantly to the observed binding kinetics, it would also be challenging to accurately predict the binding kinetics. 2.2 Details of the Simulation Protocols The covalent FEP method12 is implemented in the FEP+ program in Schrodinger Suite 2018 and used for all the relative binding free energy calculations reported in the present paper. OPLS3 force field20 was used for molecular modeling of all ligands and proteins and SPC potential is used for water. The protein systems were prepared using the Protein


10

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


Preparation Wizard21 during which protonation states were assigned assuming that pH of experimental assay is 7.0. The systems were solvated in a water box that has buffer width of 5.0 Å and 10.0 Å for complex simulation and solvent simulation. The whole systems were relaxed and equilibrated with the following relaxation protocol. First the system with restraints to its initial position was simulated by Brownie dynamics at 10 K using an NVT ensemble for 100 pico-second (ps) followed by an NPT ensemble. Second the system with restraints was simulated at 300 K using NPT ensemble, then the system without restraints was simulated at 300 K using NPT ensemble for 240 ps. Eventually the system was simulated by lambda hopping with 12 lambda windows for 10 nano-seconds (ns). Replica exchange between neighboring lambda windows were attempted every 1.2 ps. The free energy was calculated by using the Bennett acceptance ratio method (BAR)22. The sampling errors and convergence were evaluated by cycle closure analysis.23,24 The transition state energy barrier for the reaction between the ligand warhead and a suitable model of the reactive protein residue were calculated by using the B3LYP-D3 density functional in conjunction with the 6-311+G* basis set25,26.

Transition state

geometries were located using our automated transition state finding methodology, AutoTS, which only takes in structures of the reactants and products from the user and generate a transition state (TS) by performing the following steps. First, a guess at the TS, either by using predefined template or by finding a geometry corresponding to the maximum (DFT) energy along a linear interpolated path between reactants and product. Second, this guessed geometry is optimized by applying the hybrid quadratic synchronous transition (QST) and eigenvector following approach of Peng and Schlegel.27 Third, the TS structure obtained from step two is being vetted by performing the intrinsic reaction coordinate (IRC)28 calculation to ensure that the transition state found is able to connect the reactants and product along the minimum reaction pathway. Finally, if any above steps failed, a new transition state will be automatically generated and continue the workflow. A more detailed description of AutoTS algorithm is available in reference 11. The relative reaction barrier heights were corrected for solvation effects by performing single point energy calculations utilizing the SM829 implicit solvent model. Flanagan et al15 have demonstrated that this methodology is capable of reproducing relative half lives of acrylamides reacting with glutathione, a common measure of intrinsic reactivity for this reaction type, although a


11


Page 12 of 39

better treatment using range-separated functional would lead to a more stable intermediate state as pointed out by Smith et al. 30 All DFT calculations utilized the pseudospectral approximation and were performed by the Jaguar program in Schrodinger Suite 2018.31 3. Results and Discussions 3.1 BTK BTK, a nonreceptor cytoplasmic tyrosine kinase belonging to the Tec kinase family, plays an important role in B-cell proliferation, differentiation, and signaling.32 BTK is expressed in specific cells of the myeloid lineage, and is shown to contribute to immunecomplex mediated activation of the FcγR and FcεR signaling pathways.33-35 In humans, defects in the BTK gene lead to a decrease in B cell numbers and immunodeficiency. Since B cells play an essential role in regulating immune response, BTK inhibitors could be effective in treating autoimmune diseases such as rheumatoid arthritis and systemic lupus erythematous. BTK is a member of a group of eleven tyrosine kinases that contain a conserved cysteine residue adjacent to the ATP-binding site,36,37 which in turn has led to extensive efforts to design covalent inhibitors that might form a covalent linkage with this conserved cysteine residue. For many of these covalent inhibitors, the barrier height for the reaction forming the covalent linkage between the inhibitor and protein residue is relatively high, and thus such covalent inhibitors typically displayed a time-dependent inhibition profile.38-40 Further, it is well known that selective binding of BTK inhibitors is very important due to the fact that the lack of selectivity over other kinases could have potentially negative effects.41 In the present paper, we are aiming to computationally estimate the relative binding kinetics of a series of BTK irreversible covalent inhibitors using the method presented in Section 2. In the work of Kim, Maderna, et al.,13 a set of 40 compounds with an imidazoquinoxaline core were studied with the intent of advancing drug development efforts for rheumatoid arthritis. Ligands 36 and 42 (which have a fluorine or a chlorine atom at the C7 position of the fused ring respectively) showed good BTK inhibitory activity and selectivity over LYN kinase inhibition. In Figure 3, we showed the transition state structures of Compounds 36 and 42 for the reactions with the entire BKT protein and a


12

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


single cysteine residue respectively. In addition to the fluoro to chloro mutation mentioned above, Ligand 42 has an additional fluorine atom at the meta position of the terminal phenyl ring. Experimentally, the kinact / Ki value is measured as 53,200 and 11,900 M-1s-1 for Compounds 36 and 42 respectively. The reaction between the cysteine residue and the acrylamide of the ligand warhead occurs in three steps; deprotonation of cysteine, the addition of thiolate and acrylamide that forms carbanion intermediate, and the proton transfer to carbanion intermediate that form the thioether product.30 The thiolate addition step is the rate limiting step of the entire reaction. Using the protocol described in Section 2, which combines DFT calculations of the transition state barriers for the thiolate addition reaction between the ligands and a truncated cysteine residue (methane thiolate) and the FEP calculations to evaluate the contribution from protein environment to stabilize or destabilize the transition state, we 36 36 42 42 were able to calculate the ratio of binding kinetics (kinact / Ki ) / (kinact / Ki ) to be 8.2,

in agreement with the experimental ratio of 4.4. To better understand the physical factors that contribute to the experimentally observed binding kinetics, in Table 1 we decompose the calculated relative binding kinetics into individual components. If we only consider the relative transition state barrier heights for the reaction between the ligand and methane thiolate, as computed by DFT, ligand 36 is predicted to react 58.4 times faster than Ligand 42, roughly one order of magnitude larger than the experimental results. However, the FEP calculations between the two ligands in the transition state indicate that the protein environments will reduce the transition state energy barrier by 1.16 kcal/mol for Ligand 42 as compared to Ligand 36. Taking both factors into account, Ligand 36 is predicted to react 8.2 times faster than Ligand 42, which agrees with experimental result within 2 fold. The raw free energy and reaction barrier calculation results are reported in Table S1 and S4 of the supporting information. To gauge the importance of the irreversible nature of the binding process, In Table 1 we also showed the FEP results between the two ligands using the product state structures


13


Page 14 of 39

instead of the transition state structures. This treatment assumes the two ligands are reversible covalent inhibitors, and calculates their relative binding free energies. Using this protocol, Ligand 42 is calculated to bind stronger than ligand 36 by -0.37 kcal/mol and the calculated binding kinetic ratio is 0.5, which contradicts with experimental results. Therefore, we find evidence in this case that it is important to consider the irreversible nature of the binding process, and calculate the binding kinetics through modeling the transition state. In another study, Liang et al. reported a highly selective irreversible BTK inhibitor, Compound 9 (CHMFL-BTK-01) which is shown in Table 2, that displayed a high selectivity profile in KINOME scan among 468 kinases/mutants at the ligand concentration of 1μM.42 In comparison, four derivatives of Compound 9 were also made that showed lower selectivity compared to Compound 9 with measured kinact / Ki values shown in Table 2. All these five ligands react with residue Cys481 via Michael addition. These ligands (denoted as Compounds 9, 19, 30, 3, and 4) have kinact / Ki values ranging from 0.01 to 0.42 µM-1s-1. Crystal structure of Compound 9 in complex with BTK indicates residues Met477 and Ala478 form hydrogen bonds with aminopyrazinone moiety in the ligands and residues Leu542, Phe413, and Val546 interact with tert-butyl moiety of Compound 9 through hydrophobic interactions.42 Following the same protocol as the above example, we calculated the relative reaction rates of these five ligands. For these five BKT inhibitors, the perturbation region is far from the warhead, and only the FEP results in transition state structures are used to predict the relative reaction rates assuming the intrinsic reactivity is the same among these ligands. Since Compounds 3 and 4 have dramatically different chemical structures than the other three compounds (9, 19, and 30) and direct FEP calculations between ligand with dramatically different structures like this is very challenging to converge, these five compounds are separated into two groups, and relative binding kinetics is calculated between compounds within each group. In particular, the relative binding kinetics of Compound 3 is calculated in reference with Compound 4, while the relative binding


14

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


kinetics of Compounds 9 and 19 is calculated in reference with Compound 30. The irreversible covalent FEP calculated results were found to agree well with the experimental data. The error of the calculation was quantified by the difference of the calculated Δlog( kinact / Ki ) value versus the experimental value for each pair of ligands, e.g.,

k (1) / Ki(1) k (1) / Ki(1) Error  log( inact )experiment  log( inact ) prediction (2) (2) kinact / Ki(2) kinact / Ki(2)

(5)

The overall mean assigned error (MUE) for all pairs of ligands of the calculated Δlog( kinact / Ki ) as compared to experimental data for these BTK inhibitors including Compounds 36 and 42 from the previous BTK system is 0.57 log units. As a comparison, we also performed FEP simulations using the protein ligand covalent product structures (assuming reversible covalent binding and calculating their relative binding free energies), and the result has a much larger MUE of 0.84 log unit as compared to 0.57 log unit using irreversible covalent FEP, indicating the importance of correctly modeling the transition state structures during the free energy analysis. 3.2 FAAH Endocannabinoids are lipid signaling molecules that regulate a wide range of mammalian behaviors, including pain, inflammation, and cognitive or emotional state.43 It is composed of two G-coupled proteins: CB1 and CB2. Tetrahydrocannabiol (THC), the psychoactive substance in marijuana44 is able to activate CB1 and has been found to be beneficial for pain relief. However, THC has significant side effects including impaired cognition and motor function, limiting its broad clinical utility. On the other hand, inhibitors of the principal AEA(anandamine)-degrading enzyme fatty acid amide hydrolase (FAAH) have demonstrated efficacy for analgesia, anti-inflammation, anxiolysis, and antidepression without disrupting mobility, cognition or body temperature45-48, which suggests highly selective FAAH inhibitors might augment endocannabinoid signaling in vivo. In 2009, Ahn et al reported a highly selective FAAH irreversible covalent inhibitor PF-3845 with an unprecedented combination of in vivo activity and selectivity, which has shown


15


Page 16 of 39

efficacy in animal model of inflammatory pain.43 The experimental data including kinact ,

Ki and kinact / Ki are measured as 0.0033 ± 0.0002 s-1, 0.23 ± 0.03 μM, and 14,310 M-1s-1. In the present work we have studied 15 analogs of PF-3845 with experimentally measured values which were reported in reference 43. All of the ligand structures are shown in Table 3, where structural modifications are explored at the meta position of the phenyl ring in the left region. These ligands irreversibly bind to the protein target by forming a carbon-oxygen bond with serine 241 through a carbamylation reaction. By performing the DFT calculations, we found that the rate limiting step involves nucleophilic attack of the alcohol group of serine on the carbonyl carbon of the ligand accompanied by a concerted proton transfer from serine to the amine of the leaving group. Both the transition state and product structures are shown in Figure 4. In Table 3, we showed the predicted ratio of the relative reaction rates of these 15 ligands together with their experimental values. The experimentally measured values range from 17 to 12,600 M-1s-1 with different R-groups on the far-left tail of the ligand and heteroatoms mutations on the X and Y positions. The irreversible covalent FEP method is able to predict the correct trend across the 15 ligands; Ligand 7 and Ligand 15 are predicted to have the fastest and slowest reaction rate respectively, in agreement with experimental results. To quantify the accuracy of our method, we have also calculated the error by using Equation 5. The overall MUE of Δlog( kinact / Ki ) is 0.94 for all ligand pairs, with a R2 value of 0.75 between the predicted and the experimental log( kinact / Ki ) values for all the ligands as plotted in Figure 5. While a significant portion of the correlation comes from ligand 15 (bottom left point of the scatter point), the MUE of Δlog( kinact / Ki ) are roughly the same with and without ligand 15 (0.94 for all ligand pairs versus 0.81 for ligands pairs without ligand 15). 3.3 KRAS KRAS protein acts as an on/off switch in cell signaling. KRAS mutations can cause cells to proliferate out of control, and have wide-spread prevalence in human cancers.49 For


16

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


example, KRAS mutant G12C mutations are found in non-small cell lung cancer, which comprises 11%-16% of lung adenocarcinomas, as well as 1%-4% of pancreatic and colorectal adenocarcinomas, respectively. 50-53 Despite the importance of KRAS mutations, a decade of extensive efforts to develop small molecule therapeutics targeting mutated KRAS remain unsuccessful.54 In early 2018, Janes et al. reported the design and characterization of a class of S-IIP G12C inhibitors with improved potency and pharmacologic properties.55 The lead compound ARS-1620 was able to achieve sufficient in vivo covalent target occupancy to inhibit KRAS-GTP binding in tumors. In the present paper, we studied relative binding kinetics of six ligands reported in the Janes’ paper. The structures of these ligands are shown in Table 4 where the modifications explored are simple R-group additions to the fused and phenyl rings. The ligand reacts with cysteine residue via Michael addition. The structures of the protein-ligand (ligand 806) complex in transition state and in product state were shown in Figure 6. The irreversible covalent FEP protocol is able to predict the correct relative reaction rate ratio, with an overall MUE of 0.38 for the calculated Δlog( kinact / Ki ). Interestingly, for this series of ligands, the reversible covalent FEP protocol, which uses the product of protein-ligand covalent complex structure instead of the transition state structure also yielded very similar results with a MUE of 0.40 for the calculated Δlog( kinact / Ki ). As shown in Figure 6, for this particular reaction, the transition state structure is similar to the product structure, and the perturbation is relatively far from the warhead. Therefore, the interaction between the ligands and the protein is similar in the transition state versus that observed in the product. 3.4 Summary of Retrospective Studies We have studied the relative binding kinetics of four irreversible covalent inhibitor series for multiple targets of high medical interest, including BTK, FAAH, and KRAS. By applying the irreversible covalent FEP method presented in the Section 2, we were able to accurately model the relative binding kinetics among various ligands, with


17


Page 18 of 39

an overall MUE of Δlog( kinact / Ki ) of 0.79 for all investigated ligand pairs. We also observed that, for perturbations close to the ligand warhead, both the intrinsic reactivity of the warhead and the interactions between the ligand and protein in the transition state have a strong effect on the reaction rate. For example, for BTK ligands 36 and 42, the irreversible FEP method gives a small error of 0.27 log units with correct sign, where in contrast the reversible FEP method gives a much larger error of 0.92 log units with incorrect sign. On the other hand, if the perturbations are far away from the warhead and the structures of the transition state and product are very similar, FEP simulations using the transition state or the product structure are observed to both well capture the experimentally observed reaction rate. For example, for the KRAS system the irreversible and reversible covalent FEP scoring methods provide essentially identical results. 4. Conclusions The major objective of irreversible covalent ligand lead optimization is to achieve highly potent and selective binding while maintaining other properties required for safety and biological efficacy. The design of selective irreversible covalent inhibitors requires the fine-tuning of the intrinsic reactivity, as well as optimization of the interactions between the ligand and protein in the transition state associated with the formation of the covalent linkage. The most important property to optimize in irreversible covalent ligand optimization is the apparent rate of formation of the covalent complex, namely kinact / Ki In this paper, we have introduced a novel method to calculate the relative binding kinetics of irreversible covalent inhibitors. The proposed method combines quantum mechanics calculations of the transition state barrier height for the reaction between the warhead of the inhibitor and a single protein residue, and molecular mechanics-based free energy calculations to account for the effect of protein environment in the binding site on the transition state barrier. We have here applied this method to 28 ligands and 4 biological target proteins of high pharmaceutical interest. The calculated relative binding kinetics were found to agree well with experimental data, with an overall MUE of Δlog(

kinact / Ki ) of 0.79 log units for all ligand pairs. The modelled transformations include modifications near and far from the ligand warhead. We anticipate this method can


18

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


greatly help the design and optimization of irreversible covalent inhibitors in future studies. 5. Supporting Information The Supporting Information is available free of charge on the ACS Publications website at DOI: The following information is included in the SI: FEP simulation graph and free energies of BTK, FAAH, and KRAS systems, transition state geometries of ligands generated by both QM and QM/MM methods. 6. Acknowledgement The author thanks Wei Chen, Yuqing Deng, Mats Svensson, Art Bochevarov, Jeremy Rupert Greenwood, Goran Krilov, Fiona McRobb, Andrea Bortolato, and Hidenori Takahashi for the general discussions during the development of this project.


19


Page 20 of 39

Table 1. Comparison of the calculated ratioa of the kinetic rate between ligands 36 and 42 binding to BTKb from various methods versus the experimental value. The calculations were done using, FEP using the transition state (TS) structure, FEP using the product state structure, the reaction barrier (BH) differenceb computed by density functional method, and combing FEP with TS and BH difference, respectively.

Compound 36

Compound 42

Property/Method

Exp.

FEP(TS)

FEP(Prod).

DFT

DFT+FEP(TS)

Rate Ratioa

4.4

0.1

0.5

58.4

8.2

36 36 42 42 a. The rate ratio is defined as (kinact / Ki ) : (kinact / Ki )

b. The protein structure is obtained from Protein Database Bank with PDB ID of 3T9T.


20

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


Table 2 Comparison of the Predicted and Experimental relative binding kineticsa of two congeneric sets of BTK irreversible covalent ligands.

Compounds (3,4) BTK Inhibitorsb

Compounds (9,19,30)

Protein PDBc

R-groups

Compound 3

log(Pred. Ratio)d

log(Exp. Ratio)

0.49

1.07

0.0

0.0

R1= Compound 4

R1=


21

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

Compound 9

5j87

Page 22 of 39

-1.58

-0.58

-1.31

-0.92

0.0

0.0

R1= Compound 19 R2= Compound 30

R2= MUE

0.57

1 1 2 2 a. The relative kinetics is defined as Δlog( kinact / Ki ), which equals to log(kinact / Ki / kinact / Ki ) , where index 1 and 2

stands for a pair of compounds. b. Five inhibitors are placed into two groups, (3, 4) and (9, 19, 30) based on chemical similarity. c. The protein structure is obtained from Protein Database Bank with PDB ID of 5j87. d. The ratio between Compounds 3 and 4 is calculated by comparing the overall reaction rate ( kinact / Ki ) of Compound 3 to overall reaction rate for Compound 4 and the ratios among Compounds 9, 19, and 30 are calculated by comparing the overall reaction rate ( kinact / Ki ) of Compounds 9 and 19 to overall reaction rate for Compound 30. 1 1 2 2 e. The mean unsigned error is calculated by comparing log(kinact / Ki / kinact / Ki ) of all pairs of ligands in the FEP graphs


22

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


with corresponding experimental results. The FEP graphs including each pair is provided in the supporting information.


23


Page 24 of 39

Table 3. Comparison of the Predicted and Experimental relative kineticsa of 15 FAAH system irreversible covalent ligands.

Ligands (1-15) X Y Pred. Ratiob Exp. Ratio log(Pred. Ratio) log(Exp. Ratio)

Name

R

FAAH_lig1

4-MePh

C N

105.7

84.7

2.02

1.93

FAAH_lig2

4-OMePh

C N

1513.3

80.6

3.18

1.91

FAAH_lig3

4-CF3Ph

C N

4328.1

239.4

3.64

2.38

FAAH_lig4

4-CF3-2-pyr

C N

520.2

244.7

2.72

2.39

FAAH_lig5

Ph

C C

894.8

75.3

2.95

1.88

FAAH_lig6

4-MePh

C C

10448.7

272.9

4.02

2.44

FAAH_lig7

4-CF3Ph

C C

580233.8

542.9

5.76

2.73

FAAH_lig8

4-FPh

C C

2647.5

203.5

3.42

2.31

FAAH_lig9

2-pyr

C C

93.9

64.1

1.97

1.81

FAAH_lig10

4-CF3-2-pyr

C C

26540.9

741.2

4.42

2.87

FAAH_lig11

3-CF3-2-pyr

C C

424.5

15.7

2.63

1.2

FAAH_lig12

6-CF3-2-pyr

C C

648.4

25.8

2.81

1.41

26094.8

281.2

4.42

2.45

FAAH_lig13 4-OMe-2-pyr C C


24

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


FAAH_lig14 4-NH2-2-pyr C C

34

21.6

1.53

1.34

FAAH_lig15

1

1

0

0

4-MePh

N N

MUEc

0.94

1 1 2 2 a. The relative kinetics is defined as Δlog( kinact / Ki ), which equals to log(kinact / Ki / kinact / Ki ) , where index 1 and 2 stand for

a pair of compounds. b. The ratio is calculated by comparing the overall reaction rate ( kinact / Ki ) of each ligand to the overall reaction rate to ligand 15, it is defined as rate(i)/rate(15). 1 1 2 2 c. The mean unsigned error is calculated by comparing log(kinact / Ki / kinact / Ki ) of all pairs of ligands in the FEP graph with

corresponding experimental results. The FEP graph including each pair is provided in the supporting information.


25


Page 26 of 39

Table 4. Comparing Predicted and Experimental relative kineticsa of 6 KRAS covalent ligands.

Ligands (ARS-806, 869, 917, 1116, 1170, 1323) Name

R1

R2

R3

R4 Pred. Ratiob Exp. Ratio log(Pred. Ratio) log(Exp. Ratio)

ARS-806

H

Cl

H

H

1.0

1.0

0.00

0.00

ARS-869

H

Cl

OH

H

1.6

1.2

0.20

0.08

ARS-917

H

Cl

F

H

6.3

2.2

0.80

0.34

ARS-1116

F

Cl

OH

H

6.1

7.5

0.79

0.88

ARS-1170

F

Cl

H

F

10.2

11.0

1.01

1.04

Cl

F

F

11.5

35.0

1.06

1.54

ARS-1323 OH MUEc

0.38

1 1 2 2 a. The relative kinetics is defined as Δlog( kinact / Ki ), which equals to log(kinact / Ki / kinact / Ki ) , where index 1 and 2

stands for a pair of compounds. b. The ratio is calculated by comparing the overall reaction rate ( kinact / Ki ) of each ligand to the overall reaction rate of


26

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


ligand ARS-806). 1 1 2 2 c. The mean unsigned error is calculated by comparing log(kinact / Ki / kinact / Ki ) of all pairs of ligands in the FEP graph

with corresponding experimental results. The FEP graph including each pair is provided in the supporting information.


27


Page 28 of 39

(1)

Ki  koff / kon (2) Figure 1. The potential energy curve for irreversible covalent ligand binding to the target protein and the corresponding binding kinetics. kon measures the association rate of the ligand and the receptor (on-rate), koff measure the dissociation rate of the ligand and the receptor (off-rate), kinact measures the rate of covalent inactivation of the receptor and Ki is the equilibrium constant for protein ligand binding. P, L, P.L, P-L≠, and P-L stand for protein, ligand, protein-ligand complex in reactant state, protein-ligand complex in transition state, and protein-ligand complex in product state. ΔG≠ is standard Gibbs free energy of activation, ΔGi is the binding free energy, and ΔGapp≠ is the apparent free energy of activation.


28

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


Figure 2. Potential Energy Curves of Reaction between Covalent Ligand and Corresponding Residue in the Protein Binding Pocket. P+L, P.L, P-L≠, and P-L stands for the separate form of protein and ligand, protein ligand complex, transition state formed by protein and covalent ligand, product formed by protein and ligand through covalent bond, respectively. R+L, R-L≠, and R-L stand for the separate form of residue and ligand, transition state formed by residue and ligand, and product formed by residue and ligand, respectively. The overall reaction barrier E(P-L≠) – E(P+L) can be decomposed into two parts: the transition state barrier described by DFT for the reaction between the ligand and a single protein residue, and the interaction between protein environment and the ligand in transition state described by FEP.


29


Page 30 of 39

Figure 3. Thermodynamic cycle for covalent Ligands 42 and 35 binding to BTK protein target. Starting from left two pictures to the middle, Ligand 36/42 and Residue CYS442 forms transition states in solvent phase; the Gibbs activation free energy difference G1(36)  G1(42) is calculated by DFT method. From the middle two pictures to the right, Ligand 36/42 and CYS422 complex in transition state is transferred to Ligand36/42 and protein complex in transition state to account for the effect of protein environment on the transition state energy barrier; the relative effect of protein environment on the Gibbs activation free energy G 2 (36)  G 2 (42) can be calculated by FEP method.18,23,24


30

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


Figure 4. The binding poses of ligand 10 (PF-3845) in the FAAH system. Above panel represents the reaction between serine and ligand and below panel represents the protein and ligand complex in transition and product state. The key bond distances of transition state geometry are shown in the figure (unit: Å).


31


Page 32 of 39

Figure 5. Comparison of the calculated and experimental binding kinetics of 15 ligands forming covalent bonds with residue serine 241 in FAAH, the relative ratios of kinact / Ki between pairs of ligands from calculation are converted to the absolute values of

kinact / Ki for all ligands by setting the average of the prediction to be equal to the average of experimental values. The correlation (R2) is 0.75 between the predicted log( kinact / Ki ) and experimental log( kinact / Ki ).


32

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


Figure 6. The binding poses of ligand 806 in the KRAS system. Left panel represents the protein and ligand in transition state and right panel represents the protein and ligand in product state. The distances of hydrogen bonds between LYS16 and ligand, carbon-sulphur bonds are shown in the figure (unit: Å).


33


Page 34 of 39

 REFERENCES

(1) Discovery and Development of Penicillin. Am. Chem. Soc. Retrieved 30 August 2015. (2) Siddique, A.; Butt, M.; Shantsila, E.; Lip, G. Y. New antiplatelet drugs: beyond aspirin and clopidogrel. Int. J. Clin. Pract. 2009, 63, 776-789. (3) Hirsh, V. Next-Generation Covalent Irreversible Kinase Inhibitors in NSCLC: Focus on Afatinib. BioDrugs. 2015, 29, 167-183. (4) Davids, M. S.; Brown, J. R. Ibrutinib: a first in class covalent inhibitor of Bruton's tyrosine kinase. Future Oncol. 2014, 10, 957-967. (5) Barnes, T. A.; O'Kane, G. M.; Vincent, M. D.; Leighl, N. B. Third-Generation Tyrosine Kinase Inhibitors Targeting Epidermal Growth Factor Receptor Mutations in Non-Small Cell Lung Cancer. Front Oncol. 2017, 7, 113-113. (6) Proudfoot, A. T.; Wright, N. Acute paracetamol poisoning. Br Med J 1970, 3, 557-558. (7) De Cesco, S.; Kurian, J.; Dufresne, C.; Mittermaier, A. K.; Moitessier, N. Covalent inhibitors design and discovery. Eur. J. Med. Chem. 2017, 138, 96-114. (8) Kuhn, B.; Tichý, M.; Wang, L.; Robinson, S.; Martin, R. E.; Kuglstatter, A.; Benz, J.; Giroud, M.; Schirmeister, T.; Abel, R.; Diederich, F.; Hert, J. Prospective Evaluation of Free Energy Calculations for the Prioritization of Cathepsin L Inhibitors. Journal of Medicinal Chemistry 2017, 60, 2485-2497. (9) Chatterjee, P.; Botello-Smith, W. M.; Zhang, H.; Qian, L.; Alsamarah, A.; Kent, D.; Lacroix, J. J.; Baudry, M.; Luo, Y. Can Relative Binding Free Energy Predict Selectivity of Reversible Covalent Inhibitors? J. Am. Chem. Soc. 2017, 139, 17945-17952. (10) Berne, B. J.; Borkovec, M.; Straub, J. E. Classical and modern methods in reaction rate theory. J. Phys. Chem. 1988, 92, 3711-3725. (11) Jacobson, L. D.; Bochevarov, A. D.; Watson, M. A.; Hughes, T. F.; Rinaldo, D.; Ehrlich, S.; Steinbrecher, T. B.; Vaitheeswaran, S.; Philipp, D. M.; Halls, M. D.; Friesner, R. A. Automated Transition State Search and Its Application to Diverse Types of Organic Reactions. J. Chem. Theory Comput. 2017, 13, 5780-5797. (12) Kuhn, B.; Tichý, M.; Wang, L.; Robinson, S.; Martin, R. E.; Kuglstatter, A.; Benz, J.; Giroud, M.; Schirmeister, T.; Abel, R.; Diederich, F.; Hert, J. Prospective Evaluation of Free Energy Calculations for the Prioritization of Cathepsin L Inhibitors. J. Med. Chem. 2017, 60, 2485-2497. (13) Kim, K.-H.; Maderna, A.; Schnute, M. E.; Hegen, M.; Mohan, S.; Miyashiro, J.; Lin, L.; Li, E.; Keegan, S.; Lussier, J.; Wrocklage, C.; Nickerson-Nutter, C. L.; Wittwer, A. J.; Soutter, H.; Caspers, N.; Han, S.; Kurumbail, R.; Dunussi-Joannopoulos, K.; Douhan, J.; Wissner, A. Imidazo[1,5-a]quinoxalines as irreversible BTK inhibitors for the treatment of rheumatoid arthritis. Bioorg. Med. Chem. Lett. 2011, 21, 6258-6263. (14) Lonsdale, R.; Ward, R. A. Structure-based design of targeted covalent inhibitors. Chem. Soc. Rev. 2018, 47, 3816-3830. (15) Flanagan, M. E.; Abramite, J. A.; Anderson, D. P.; Aulabaugh, A.; Dahal, U. P.; Gilbert, A. M.; Li, C.; Montgomery, J.; Oppenheimer, S. R.; Ryder, T.; Schuff, B. P.; Uccello, D. P.; Walker, G. S.; Wu, Y.; Brown, M. F.; Chen, J. M.; Hayward, M. M.; Noe, M. C.; Obach, R. S.; Philippe, L.; Shanmugasundaram, V.; Shapiro, M. J.; Starr, J.; Stroh, J.; Che, Y. Chemical and Computational Methods for the Characterization of Covalent Reactive Groups for the Prospective Design of Irreversible Inhibitors. J. Med. Chem. 2014, 57, 10072-10079.


34

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


(16) Abel, R.; Wang, L.; Harder, E. D.; Berne, B. J.; Friesner, R. A. Advancing Drug Discovery through Enhanced Free Energy Calculations. Acc. Chem. Res. 2017, 50, 1625-1632. (17) Wang, L.; Deng, Y.; Wu, Y.; Kim, B.; LeBard, D. N.; Wandschneider, D.; Beachy, M.; Friesner, R. A.; Abel, R. Accurate Modeling of Scaffold Hopping Transformations in Drug Discovery. J. Chem. Theory Comput. 2017, 13, 42-54. (18) Wang, L.; Wu, Y.; Deng, Y.; Kim, B.; Pierce, L.; Krilov, G.; Lupyan, D.; Robinson, S.; Dahlgren, M. K.; Greenwood, J.; Romero, D. L.; Masse, C.; Knight, J. L.; Steinbrecher, T.; Beuming, T.; Damm, W.; Harder, E.; Sherman, W.; Brewer, M.; Wester, R.; Murcko, M.; Frye, L.; Farid, R.; Lin, T.; Mobley, D. L.; Jorgensen, W. L.; Berne, B. J.; Friesner, R. A.; Abel, R. 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. (19) Rydberg, P.; Olsen, L.; Norrby, P. O.; Ryde, U. General Transition-State Force Field for Cytochrome P450 Hydroxylation. J Chem Theory Comput 2007, 3, 1765-1773. (20) Harder, E.; Damm, W.; Maple, J.; Wu, C.; Reboul, M.; Xiang, J. Y.; Wang, L.; Lupyan, D.; Dahlgren, M. K.; Knight, J. L.; Kaus, J. W.; Cerutti, D. S.; Krilov, G.; Jorgensen, W. L.; Abel, R.; Friesner, R. A. OPLS3: A Force Field Providing Broad Coverage of Drug-like Small Molecules and Proteins. J. Chem. Theory Comput. 2016, 12, 281-296. (21) In Schrodinger Suite 2019 Protein Preparation Wizard. Schrodinger L.L.C: New York, NY, 2019. (22) Bennett, C. H. Efficient estimation of free energy differences from Monte Carlo data. J. Comput. Phys. 1976, 22, 245-268. (23) Wang, L.; Deng, Y.; Knight, J. L.; Wu, Y.; Kim, B.; Sherman, W.; Shelley, J. C.; Lin, T.; Abel, R. Modeling Local Structural Rearrangements Using FEP/REST: Application to Relative Binding Affinity Predictions of CDK2 Inhibitors. J. Chem. Theory Comput. 2013, 9, 1282-1293. (24) Lingle Wang, T. L., Robert Abel. Cycle Closure Estimation of Relative Binding Affinities and Errors. (25) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (26) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623-11627. (27) Peng, C.; Bernhard Schlegel, H. Combining Synchronous Transit and Quasi-Newton Methods to Find Transition States. abrégé en Isr. J. Chem. 1993, 33, 449-454. (28) Gonzalez, C.; Schlegel, H. B. An improved algorithm for reaction path following. J. Chem. Phys. 1989, 90, 2154-2161. (29) Cramer, C. J.; Truhlar, D. G. A Universal Approach to Solvation Modeling. Acc. Chem. Res. 2008, 41, 760-768. (30) Smith, J. M.; Jami Alahmadi, Y.; Rowley, C. N. Range-Separated DFT Functionals are Necessary to Model Thio-Michael Additions. J. Chem. Theory Comput. 2013, 9, 4860-4865. (31) Bochevarov, A. D.; Harder, E.; Hughes, T. F.; Greenwood, J. R.; Braden, D. A.; Philipp, D. M.; Rinaldo, D.; Halls, M. D.; Zhang, J.; Friesner, R. A. Jaguar: A highperformance quantum chemistry software program with strengths in life and materials sciences. Int. J. Quantum Chem. 2013, 113, 2110-2142.


35


Page 36 of 39

(32) Chang, B. Y.; Huang, M. M.; Francesco, M.; Chen, J.; Sokolove, J.; Magadala, P.; Robinson, W. H.; Buggy, J. J. The Bruton tyrosine kinase inhibitor PCI-32765 ameliorates autoimmune arthritis by inhibition of multiple effector cells. Arthritis Res. Ther. 2011, 13, R115. (33) Hata, D.; Kawakami, Y.; Inagaki, N.; Lantz, C. S.; Kitamura, T.; Khan, W. N.; MaedaYamamoto, M.; Miura, T.; Han, W.; Hartman, S. E.; Yao, L.; Nagai, H.; Goldfeld, A. E.; Alt, F. W.; Galli, S. J.; Witte, O. N.; Kawakami, T. Involvement of Bruton's tyrosine kinase in FcepsilonRI-dependent mast cell degranulation and cytokine production. J Exp Med 1998, 187, 1235-1247. (34) Jongstra-Bilen, J.; Puig Cano, A.; Hasija, M.; Xiao, H.; Smith, C. I.; Cybulsky, M. I. Dual functions of Bruton's tyrosine kinase and Tec kinase during Fcgamma receptor-induced signaling and phagocytosis. J. Immunol. 2008, 181, 288-298. (35) Wang, X. S.; Lau, H. Y. Histamine release from human buffy coat-derived mast cells. Int Immunopharmacol. 2007, 7, 541-546. (36) Cohen, M. S.; Zhang, C.; Shokat, K. M.; Taunton, J. Structural Bioinformatics-Based Design of Selective, Irreversible Kinase Inhibitors. Science 2005, 308, 1318-1321. (37) Leproult, E.; Barluenga, S.; Moras, D.; Wurtz, J.-M.; Winssinger, N. Cysteine Mapping in Conformationally Distinct Kinase Nucleotide Binding Sites: Application to the Design of Selective Covalent Inhibitors. J. Med. Chem. 2011, 54, 1347-1355. (38) Singh, J.; Dobrusin, E. M.; Fry, D. W.; Haske, T.; Whitty, A.; McNamara, D. J. Structure-based design of a potent, selective, and irreversible inhibitor of the catalytic domain of the erbB receptor subfamily of protein tyrosine kinases. J Med Chem 1997, 40, 1130-1135. (39) Gierse, J. K.; Koboldt, C. M.; Walker, M. C.; Seibert, K.; Isakson, P. C. Kinetic basis for selective inhibition of cyclo-oxygenases. Biochem J 1999, 339 ( Pt 3), 607-614. (40) Wang, Z.; Watt, W.; Brooks, N. A.; Harris, M. S.; Urban, J.; Boatman, D.; McMillan, M.; Kahn, M.; Heinrikson, R. L.; Finzel, B. C.; Wittwer, A. J.; Blinn, J.; Kamtekar, S.; Tomasselli, A. G. Kinetic and structural characterization of caspase-3 and caspase-8 inhibition by a novel class of irreversible inhibitors. Biochim Biophys Acta 2010, 1804, 1817-1831. (41) Patel, V.; Balakrishnan, K.; Bibikova, E.; Ayres, M.; Keating, M. J.; Wierda, W. G.; Gandhi, V. Comparison of Acalabrutinib, A Selective Bruton Tyrosine Kinase Inhibitor, with Ibrutinib in Chronic Lymphocytic Leukemia Cells. Clin. Cancer Res. 2017, 23, 3734-3743. (42) Liang, Q.; Chen, Y.; Yu, K.; Chen, C.; Zhang, S.; Wang, A.; Wang, W.; Wu, H.; Liu, X.; Wang, B.; Wang, L.; Hu, Z.; Wang, W.; Ren, T.; Zhang, S.; Liu, Q.; Yun, C.-H.; Liu, J. Discovery of N-(3-(5-((3-acrylamido-4-(morpholine-4-carbonyl)phenyl)amino)-1-methyl-6oxo-1,6-dihydropyridin-3-yl)-2-methylphenyl)-4-(tert-butyl)benzamide (CHMFL-BTK-01) as a highly selective irreversible Bruton's tyrosine kinase (BTK) inhibitor. Eur. J. Med. Chem. 2017, 131, 107-125. (43) Ahn, K.; Johnson, D. S.; Mileni, M.; Beidler, D.; Long, J. Z.; McKinney, M. K.; Weerapana, E.; Sadagopan, N.; Liimatta, M.; Smith, S. E.; Lazerwith, S.; Stiff, C.; Kamtekar, S.; Bhattacharya, K.; Zhang, Y.; Swaney, S.; Van Becelaere, K.; Stevens, R. C.; Cravatt, B. F. Discovery and Characterization of a Highly Selective FAAH Inhibitor that Reduces Inflammatory Pain. Chem. Biol. 2009, 16, 411-420. (44) Mechoulam, R.: The Pharmacohistory of Cannabis sativa; CRS Press: FL, 1986. (45) Boger, D. L.; Sato, H.; Lerner, A. E.; Hedrick, M. P.; Fecik, R. A.; Miyauchi, H.; Wilkie, G. D.; Austin, B. J.; Patricelli, M. P.; Cravatt, B. F. Exceptionally potent inhibitors of fatty acid amide hydrolase: The enzyme responsible for degradation of endogenous oleamide


36

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


and anandamide. Proc. Natl. Acad. Sci. 2000, 97, 5044-5049. (46) Cravatt, B. F.; Demarest, K.; Patricelli, M. P.; Bracey, M. H.; Giang, D. K.; Martin, B. R.; Lichtman, A. H. Supersensitivity to Anandamide and Enhanced Endogenous Cannabinoid Signaling in Mice Lacking Fatty Acid Amide Hydrolase. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 9371-9376. (47) Kathuria, S.; Gaetani, S.; Fegley, D.; Valiño, F.; Duranti, A.; Tontini, A.; Mor, M.; Tarzia, G.; Rana, G. L.; Calignano, A.; Giustino, A.; Tattoli, M.; Palmery, M.; Cuomo, V.; Piomelli, D. Modulation of anxiety through blockade of anandamide hydrolysis. Nat. Med. 2002, 9, 76. (48) Lichtman, A. H.; Leung, D.; Shelton, C. C.; Saghatelian, A.; Hardouin, C.; Boger, D. L.; Cravatt, B. F. Reversible Inhibitors of Fatty Acid Amide Hydrolase That Promote Analgesia: Evidence for an Unprecedented Combination of Potency and Selectivity. J. Pharmacol. Exp. Ther. 2004, 311, 441. (49) Bos, J. L. ras oncogenes in human cancer: a review. Cancer Res. 1989, 49, 4682-4689. (50) Bailey, P.; Chang, D. K.; Nones, K.; Johns, A. L.; Patch, A.-M.; Gingras, M.-C.; Miller, D. K.; Christ, A. N.; Bruxner, T. J. C.; Quinn, M. C.; Nourse, C.; Murtaugh, L. C.; Harliwong, I.; Idrisoglu, S.; Manning, S.; Nourbakhsh, E.; Wani, S.; Fink, L.; Holmes, O.; Chin, V.; Anderson, M. J.; Kazakoff, S.; Leonard, C.; Newell, F.; Waddell, N.; Wood, S.; Xu, Q.; Wilson, P. J.; Cloonan, N.; Kassahn, K. S.; Taylor, D.; Quek, K.; Robertson, A.; Pantano, L.; Mincarelli, L.; Sanchez, L. N.; Evers, L.; Wu, J.; Pinese, M.; Cowley, M. J.; Jones, M. D.; Colvin, E. K.; Nagrial, A. M.; Humphrey, E. S.; Chantrill, L. A.; Mawson, A.; Humphris, J.; Chou, A.; Pajic, M.; Scarlett, C. J.; Pinho, A. V.; Giry-Laterriere, M.; Rooman, I.; Samra, J. S.; Kench, J. G.; Lovell, J. A.; Merrett, N. D.; Toon, C. W.; Epari, K.; Nguyen, N. Q.; Barbour, A.; Zeps, N.; Moran-Jones, K.; Jamieson, N. B.; Graham, J. S.; Duthie, F.; Oien, K.; Hair, J.; Grützmann, R.; Maitra, A.; Iacobuzio-Donahue, C. A.; Wolfgang, C. L.; Morgan, R. A.; Lawlor, R. T.; Corbo, V.; Bassi, C.; Rusev, B.; Capelli, P.; Salvia, R.; Tortora, G.; Mukhopadhyay, D.; Petersen, G. M.; Australian Pancreatic Cancer Genome, I.; Munzy, D. M.; Fisher, W. E.; Karim, S. A.; Eshleman, J. R.; Hruban, R. H.; Pilarsky, C.; Morton, J. P.; Sansom, O. J.; Scarpa, A.; Musgrove, E. A.; Bailey, U.-M. H.; Hofmann, O.; Sutherland, R. L.; Wheeler, D. A.; Gill, A. J.; Gibbs, R. A.; Pearson, J. V. Genomic analyses identify molecular subtypes of pancreatic cancer. Nature 2016, 531, 47. (51) Campbell, J. D.; Alexandrov, A.; Kim, J.; Wala, J.; Berger, A. H.; Pedamallu, C. S.; Shukla, S. A.; Guo, G.; Brooks, A. N.; Murray, B. A.; Imielinski, M.; Hu, X.; Ling, S.; Akbani, R.; Rosenberg, M.; Cibulskis, C.; Ramachandran, A.; Collisson, E. A.; Kwiatkowski, D. J.; Lawrence, M. S.; Weinstein, J. N.; Verhaak, R. G. W.; Wu, C. J.; Hammerman, P. S.; Cherniack, A. D.; Getz, G.; Cancer Genome Atlas Research, N.; Zenklusen, J. C.; Zhang, J.; Felau, I.; Demchok, J. A.; Yang, L.; Wang, Z.; Ferguson, M. L.; Tarnuzzer, R.; Hutter, C. M.; Sofia, H. J.; Pihl, T.; Wan, Y.; Chudamani, S.; Liu, J.; Sun, C.; Naresh, R.; Lolla, L.; Wu, Y.; Creighton, C. J.; Rathmell, W. K.; Auman, J. T.; Balu, S.; Bodenheimer, T.; Hayes, D. N.; Hoadley, K. A.; Hoyle, A. P.; Jones, C. D.; Jefferys, S. R.; Meng, S.; Mieczkowski, P. A.; Mose, L. E.; Perou, C. M.; Roach, J.; Shi, Y.; Simons, J. V.; Skelly, T.; Soloway, M. G.; Tan, D.; Wu, J.; Veluvolu, U.; Parker, J. S.; Wilkerson, M. D.; Boice, L.; Huang, M.; Thorne, L. B.; Getz, G.; Noble, M. S.; Zhang, H.; Heiman, D. I.; Cho, J.; Gehlenborg, N.; Saksena, G.; Voet, D.; Lin, P.; Frazer, S.; Kim, J.; Lawrence, M. S.; Chin, L.; Tsao, M.-S.; Allison, F.; Chadwick, D.; Muley, T.; Meister, M.; Dienemann, H.; Kucherlapati, R.; Park, P.; Bowen, J.; GastierFoster, J. M.; Gerken, M.; Leraas, K. M.; Lichtenberg, T. M.; Ramirez, N. C.; Wise, L. Distinct


37


Page 38 of 39

patterns of somatic genome alterations in lung adenocarcinomas and squamous cell carcinomas. Nature Genet. 2016, 48, 607. (52) Giannakis, M.; Mu, Xinmeng J.; Shukla, Sachet A.; Qian, Zhi R.; Cohen, O.; Nishihara, R.; Bahl, S.; Cao, Y.; Amin-Mansour, A.; Yamauchi, M.; Sukawa, Y.; Stewart, C.; Rosenberg, M.; Mima, K.; Inamura, K.; Nosho, K.; Nowak, Jonathan A.; Lawrence, Michael S.; Giovannucci, Edward L.; Chan, Andrew T.; Ng, K.; Meyerhardt, Jeffrey A.; Van Allen, Eliezer M.; Getz, G.; Gabriel, Stacey B.; Lander, Eric S.; Wu, Catherine J.; Fuchs, Charles S.; Ogino, S.; Garraway, Levi A. Genomic Correlates of Immune-Cell Infiltrates in Colorectal Carcinoma. Cell Rep. 2016, 15, 857-865. (53) Jordan, E. J.; Kim, H. R.; Arcila, M. E.; Barron, D.; Chakravarty, D.; Gao, J.; Chang, M. T.; Ni, A.; Kundra, R.; Jonsson, P.; Jayakumaran, G.; Gao, S. P.; Johnsen, H. C.; Hanrahan, A. J.; Zehir, A.; Rekhtman, N.; Ginsberg, M. S.; Li, B. T.; Yu, H. A.; Paik, P. K.; Drilon, A.; Hellmann, M. D.; Reales, D. N.; Benayed, R.; Rusch, V. W.; Kris, M. G.; Chaft, J. E.; Baselga, J.; Taylor, B. S.; Schultz, N.; Rudin, C. M.; Hyman, D. M.; Berger, M. F.; Solit, D. B.; Ladanyi, M.; Riely, G. J. Prospective Comprehensive Molecular Characterization of Lung Adenocarcinomas for Efficient Patient Matching to Approved and Emerging Therapies. CANCER DISCOV. 2017, 7, 596. (54) Cox, A. D.; Fesik, S. W.; Kimmelman, A. C.; Luo, J.; Der, C. J. Drugging the undruggable RAS: Mission Possible? Nat. Rev. Drug Discov. 2014, 13, 828. (55) Janes, M. R.; Zhang, J.; Li, L. S.; Hansen, R.; Peters, U.; Guo, X.; Chen, Y.; Babbar, A.; Firdaus, S. J.; Darjania, L.; Feng, J.; Chen, J. H.; Li, S.; Li, S.; Long, Y. O.; Thach, C.; Liu, Y.; Zarieh, A.; Ely, T.; Kucharski, J. M.; Kessler, L. V.; Wu, T.; Yu, K.; Wang, Y.; Yao, Y.; Deng, X.; Zarrinkar, P. P.; Brehmer, D.; Dhanak, D.; Lorenzi, M. V.; Hu-Lowe, D.; Patricelli, M. P.; Ren, P.; Liu, Y. Targeting KRAS Mutant Cancers with a Covalent G12C-Specific Inhibitor. Cell 2018, 172, 578-589 e517.


38

Page Journal 39 ofof39 Chemical Information and Modeling KRAS Ligands

1 2 3 4 5 6

FEP + AutoTS

LSY16

ACS Paragon Plus Environment CYS12

Towards Atomistic Modelling of Irreversible Covalent Inhibitor Binding

Recommend Documents