Molecular Dynamics Analysis of Binding of Kinase Inhibitors to WT

Subscriber access provided by University of Massachusetts Amherst Libraries

Article

Molecular dynamics analysis of binding of kinase inhibitors to WT EGFR and the T790M mutant Jiyong Park, Joseph J. McDonald, Russell C. Petter, and K. N. Houk J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.5b01221 • Publication Date (Web): 24 Mar 2016 Downloaded from http://pubs.acs.org on March 25, 2016

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

Journal of Chemical Theory and Computation 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 44

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 Theory and Computation

Molecular dynamics analysis of binding of kinase inhibitors to WT EGFR and the T790M mutant Jiyong Park,† Joseph J. McDonald,‡ Russell C. Petter,∗,‡,¶ and K. N. Houk∗,† Department of Chemistry and Biochemistry, University of California, Los Angeles, CA 90095, USA, and Celgene Avilomics Research, Bedford, MA 01730, USA E-mail: [email protected]; [email protected]

Abstract Epidermal growth factor receptor (EGFR) inhibitors interrupt EGFR-dependent cellular signaling pathways that lead to accelerated tumor growth and proliferation. Mutation of a threonine in the inhibitor binding pocket, known as the ‘gatekeeper,’ to methionine (T790M) confers acquired resistance to several EGFR-selective inhibitors. We studied interactions between EGFR inhibitors and the gatekeeper residues of the target protein. Thermodynamic integration (TI) with Amber14 indicates that the binding energies of gefitinib and AEE788 to the active state of the T790M mutant EGFR is 3 kcal/mol higher than to the WT, whereas ATP binding energy to the mutant is similar to the WT. Using metadynamics MD simulations with NAMD v2.9, the conformational equilibrium between the inactive resting state and the catalytically competent activate state was determined for the WT EGFR. When combined with the results obtained by Sutto and Gervasio, our simulations showed that the T790M point mutation lowers the free energy of the active state by 5 kcal/mol relative to the ∗

To whom correspondence should be addressed UCLA ‡ Celgene ¶ Current address: Arrakis Therapeutics, Boston, MA, USA †

1 ACS Paragon Plus Environment


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

inactive state of the enzyme. Relative to the WT, the T790M mutant binds gefitinib more strongly. The T790M mutation is nevertheless resistant due to its increased binding of ATP. By contrast, the binding of AEE788 to the active state causes a conformational change in the αC-helix adjacent to the inhibitor binding pocket, that results in a 2 kcal/mol energy penalty. The energy penalty explains why the binding of AEE788 to the T790M mutant is unfavorable relative to binding to WT EGFR. These results establish the role of the gatekeeper mutation on inhibitor selectivity. Additional molecular dynamics (MD) simulations, TI, and metadynamics MD simulations reveal the origins of the changes in binding energy of WT and mutants.

1

Introduction

Epidermal growth factor receptor (EGFR) kinase is a cell surface receptor protein complex that mediates extracellular growth signals. 1 Upon binding of epidermal growth factor (EGF) signaling molecules, signal transduction is mediated by dimerization and transphosphorylation of intracellular domains and subsequent activation of the kinase domain. The activated kinase domain phosphorylates downstream signaling molecules using ATP. 2 Elevated activities of EGFR are associated with the onset and the severity of cancers including non-small-cell lung cancer (NSCLC), and glioblastoma. 3 A somatic mutation (L858R) is known to stabilize the enzymatically active state of EGFR, which is associated with onset of the cancers exploiting EGFR-dependent signaling pathways. EGFR is one of prime targets for cancer inhibition and treatment. 4 There are three EGFR-selective small molecular inhibitors approved by FDA: gefitinib (brand name: Iressa), erlotinib (brand name: Tarceva) and lapatinib (brand name: Tykerb) (Fig. 1). 5,6 The therapeutic benefits of the small molecular inhibitors of EGFR kinase have been hampered by acquired drug resistance. 3,7 After 9 to 13 months of treatment with EGFR-selective inhibitors, patients develop resistance. 8 The most prevalent mechanism of the acquired drugresistance is a somatic mutation, T790M, which generally follows a cancer initiation mutation 2 ACS Paragon Plus Environment

Page 2 of 44

Page 3 of 44

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


Figure 1: EGFR kinase domain and selective inhibitors. Chemical structures of EGFR inhibitors. (L858R). 8–11 The threonine, T790, in the active site of EGFR is known as the gatekeeper residue, that significantly influences the binding of EGFR inhibitors. The T790M mutation was originally identified from patients treated with gefitinib and later found in other patients receiving erlotinib treatment. 12 Clinical studies demonstrated that the gatekeeper mutation (T790M) resulted in a dramatic reduction in therapeutic efficacy of EGFR-selective inhibitors. 13 Recently novel inhibitors are designed to target inhibitor resistant mutants including afatinib 14 and neratinib, 15 both of which have advanced to phase III clinical trials. Two competing mechanisms have been proposed to explain the origin of the acquired drug resistance induced by the T790M somatic mutation. Kobayashi et al. suggested that steric hindrance to binding occurs due to increased molecular size of methionine over threonine. 8 However gefitinib and AEE788 are cocrystalized with the T790M mutant, 16,17 which undermines the steric hindrance mechanism. Later Yun et al. demonstrated the gatekeeper mutation (T790M) increases sensitivity to the natural substrate (ATP), while the affinities of EGFR inhibitors remains in the low nanomolar range. This is supported by biochemical evidences: a double mutant (L858R/T790M) reduces the equilibrium ATP binding constant (Km ) by 18 fold relative to a single mutant (L858R), and decreases the binding constant of EGFR inhibitors by 5–18 fold. The elevated ATP sensitivity amplifies the difference of inhibitor binding affinity at cellular concentration of ATP. This explains the reduced efficacy



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

of the inhibitors. 17

Figure 2: Conformational states of EGFR kinase domain. (a) A rendition of the active state of EGFR kinase domain. Two of structural elements are highlighted in purple (αC helix) and in yellow (activation loop). (b) Substrate binding site of EGFR with bound gefitinib (PDB Code: 2ITO). (c) Src/Cdk-like inactive state of EGFR kinase domain (PDB Code: 2GS7). The molecules are rendered using PyMOL software. 18

The increased affinity of ATP to the T790M mutant is explained by the conformational equilibrium shift. The kinase domain of EGFR adopts distinct conformational states, including a catalytically active state and and an inactive (Src/Cdk-like) state (Fig. 2). 19,20 Oncogenic mutations of EGFR (T790M, L858R, and L858R/T790M) are known to stabilize the active state of EGFR, thereby increasing the binding of ATP. 10,21,22 As an inhibitor binds to a specific conformational state of the kinase domain, the conformational equilibrium shift influences the binding affinity of the inhibitor. 23 How the T790M mutation influences the binding of the EGFR selective inhibitors remains elusive. Yun et al. also discovered that the T790M single mutation reduces Kd of gefitinib by 8 fold, which contradicts the notion that the mutation hampers the action of gefitinib. 17 The increased sensitivity to ATP could account for the observed inhibitor resistance caused by the T790M mutation. Molecular dynamics (MD) simulations have been undertaken to provide a detailed understanding of the kinase binding domain of EGFR. Previous long-time (> 1 µs) MD simulations 4 ACS Paragon Plus Environment

Page 4 of 44

Page 5 of 44

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


by the Shaw group provided information on the structural transition from active to inactive states of EGFR kinase. 24 After 1 µs of MD simulation, the structural transformation from the active state to the Src/Cdk-like inactive state was observed; a portion of the active state conformation becomes disordered and then refolds to become a Src/Cdk-like inactive state. Sutto and Gervasio carried out a series of metadynamics MD simulations and explained how point and dual mutations including T790M, L858R, and T790M/L858R modulate the conformational equilibrium between the active and the Src/Cdk-like inactive conformations of the EGFR kinase domain. 25 All three mutations are predicted to shift the conformational equilibrium toward the active state of the EGFR kinase domain. The finding is in agreement with the clinical evidence that those mutations overactivate the enzyme in human tumors, a phenomenon known as constitutive activation. 17,26–28 Although these simulations provided quantitative information on the pathway and the energetics of the kinase domain, detailed information on the origins of the binding affinity changes associated with the gatekeeper mutation were not elucidated. We have studied how the gatekeeper mutation (T790M) in the active site of EGFR kinase domain alters binding of the inhibitors. The free energy minimum state of apo EGFR is the inactive state, 25 but the binding of inhibitors occurs to the active state. Therefore the equilibrium of the two conformational states modulates the binding free energy of the inhibitors. We have computed the free energy difference between the native inactive state of the enzyme and the active conformation to which the ligands bind. The metadynamics method was used for the computation of the free energy difference of the two conformational states. 29 Using the thermodynamic integration (TI) 30 and the metadynamics method, we have computed the changes of the binding energies of inhibitors for the active conformation of EGFR upon the gatekeeper mutation (T790M). The net effect of both the binding energy differences to the active conformation and the constitutive activation explains quantitatively the observed binding affinity profile of gefitinib and AEE788. In addition, we have investigated why the gatekeeper residues change conformations upon binding of inhibitors. A combination of TI



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

Page 6 of 44

and metadynamics simulations show that the aromatic group of the inhibitor interacts with the methionine gatekeeper and causes unfavorable steric interactions. Finally we demonstrate the T790M mutant reduces the conformational flexibility of αC-helix, which explains increased stability of the active state of the T790M mutant. Since acquired drug resistance induced by the gatekeeper mutation is widespread across receptor kinases associated with cancers, 31 the detailed understanding of binding energy changes caused by mutations will aid further development of therapeutics that interrupt receptor-kinase-dependent cellular signaling pathways.

2

Results and Discussion

2.1

The gatekeeper mutation (T790M) influences the binding of gefitinib and AEE788 to the active state of EGFR

We first questioned whether the T790M point mutation induces an energetic penalty to the binding of EGFR inhibitors. We studied the interaction of a reversible inhibitors, gefitinib and AEE788 32 with EGFR, since both clinical activity 33 and equilibrium binding affinity 17 have been reported previously (Table 1). ATP and two other inhibitors were explored as well. High resolution crystal structures of inhibitor/substrate-bound cocrystals indicate that the small molecules interact with the active state of wild-type (WT) and L858R and T790M mutant (Table S1). Thus we considered the change in inhibitor/substrate binding energy to the active state upon mutation of WT gatekeeper residue (T790) to methionine of EGFR. TI simulations were performed on both the protein-drug complex and on the free protein. Using a thermodynamic cycle depicted in Fig. S1, we computed the difference in binding free energy to the active state upon mutation of the gatekeeper residue;

Active Active ∆∆GActive T 790M = ∆Gbind,M 790 − ∆Gbind,T 790


(1)

Page 7 of 44

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


Table 1: Experimental equilibrium dissociation constant (Kd ) and half maximum inhibitory concentration (IC50 ) of gefitinib and AEE788. Gefitinib Mutation Kd 17 IC50 33 WT 35.3 3.0 4.6 N/A T790M L858R 2.4 0.8

AEE788 Kd IC50 34,35 5.3 2 27.6 >2000 1.1 6 17

The unit of Kd and IC50 is nM/L.

Active , where ∆GActive bind,M 790 and ∆Gbind,T 790 are the binding free energies to the mutated gatekeeper

residue (M790) and to the WT gatekeeper residue (T790), in their active state. The technique was shown to reproduce the mutation dependent binding free energy difference of a small molecular inhibitor of p38α MAP kinase. 36 Table 2: Computed relative binding energy of EGFR inhibitors to the active conformation of EGFR kinase. Energy unit is kcal/mol. Compound ∆∆GActive T 790M Gefitinib 3.9 ± 1.0 Gefitinib 2.3 ± 1.0 Gefitinib 2.6 ± 0.7 2.6 ± 0.7 AEE788 3.0 ± 0.7 AEE788 0.2 ± 0.6 ATP WZ4002 −1.3 ± 0.7 a

PDB code 2ITO 3UG2 2ITZ a 2J6M 2JIU 2EB3 3IKA

The L858R mutation is present.

Table 2 summarizes the results of TI simulations starting from different crystal structures. The T790M mutation reduces binding affinity of gefitinib and of AEE788 to the active conformation about 3 kcal/mol. The results for gefitinib involves TI simulations using three active state structures of the WT and the L858R mutant EGFR. AEE788 is an EGFR-selective inhibitor which suffers from reduced efficacy upon introduction of the T790M mutation. 17 Two crystal structures was used for the simulations. The phase I clinical trial of AEE788 was terminated due to the toxicity and the minimal therapeutic benefit. 34,37 Both gefitinib and AEE788 have a bulky substituent interacting with the gatekeeper. The T790M mutation 7 ACS Paragon Plus Environment


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

reduces the binding affinity of AEE788 by 3 kcal/mol. The orientation of the αC-helix in the crystal structure of the T790M bound to AEE788 differs from that of the WT as shown in Fig. 5. We discuss the influence of the distinct orientation of αC-helix in Sec. 2.4. We also computed ∆∆GActive T 790M of the natural substrate (ATP). The binding affinity of ATP was found to be affected weakly by T790M mutation in a previous experiment. 17 The computed ∆∆GActive T 790M of ATP is only 0.2 ± 0.6 kcal/mol, which indicates that ATP does not experience a change in steric effects upon T790M point mutation. WZ4002 is an EGFR selective inhibitor targeting cancers which express the T790M mutation. 38 The inhibitor suppresses inhibitor resistant cancer cell lines in vitro. 39 The inhibitor shows improved binding affinity to M790 over T790 in the active state (∆∆GActive T 790M = −1.3 ± 0.7 kcal/mol). ATP and WZ4002 do not have bulky groups interacting with gatekeeper residues (Fig. S8), and the gatekeeper mutation influences the binding affinity weakly. The TI simulations demonstrate the methionine gatekeeper modulates the binding affinity of EGFR-selective inhibitors. The T790M point mutation contributes to the binding free energy differences of the inhibitors, disfavoring gefitinib and AEE788 binds to the active conformation, having little effect on ATP and favoring WZ4002 binding. The binding to the active conformation does not explain the differences in binding affinities observed experimentally (Table 1), but we show below that the conformational equilibrium shift caused by the T790M mutation also influences the overall binding free energy of the EGFR inhibitors.

2.2

The constitutive activation mutations modulate the conformational equilibrium of apo EGFR

We have demonstrated that the gatekeeper mutation causes steric hindrance to certain EGFR inhibitors. In what follows we discuss how the gatekeeper mutation (T790M) causes a conformational equilibrium shift toward the active state of EGFR, a phenomenon known as ‘constitutive activation.’ 40 Structural studies show that the resting state of EGFR is the enzymatically inactive 8 ACS Paragon Plus Environment

Page 8 of 44

Page 9 of 44

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


state, whereas binding of substrates or certain mutations activate the enzyme. 41,42 Constitutive activation mutations such as L858R, T790M, and L858R/T790M overactivate EGFR by shifting the conformational equilibrium of the catalytic site toward the active conformation. Some forms of constitutive activation increase the susceptibilities of cancer cells to gefitinib. 43,44 As inhibitors targeting EGFR recognize the active state of the protein, the conformational equilibrium shift also changes the binding affinities of the inhibitors. We first studied the conformational equilibrium of WT EGFR. The free energy difference between the active state and the Src/Cdk-like inactive state is defined as:

∆GI→A = ∆GActive − ∆GInactive .

(2)

, where ∆GInactive and ∆GActive are the conformational free energies of the inactive and the active states, respectively, shown earlier in Fig. 2. Fig. 3 summarizes metadynamics simulations of WT EGFR. The Gibbs free energy is plotted as a function of the root-mean-squared (RMS) deviation from the crystal structure of the active state (RM SDActive ) and the inactive state (RM SDInactive ). The active state is predicted to be higher in energy than the Src/Cdk-like inactive state, as known experimentally. The global free energy minimum corresponds to the Src/Cdk-like inactive state (Fig. 3: I), and the active state (Fig. 3: A) is 5 kcal/mol higher in energy. The free energy difference between the active state and the Src/Cdk-like inactive state converged: the difference in ∆GI→A at 2.4 µs and 2.7 µs is lower than 1.0 kcal/mol, as shown in Fig. S6. The local minimum corresponding to the active state in Fig. 3 overlaps with the region sampled from a 200 ns unbiased MD simulation of WT EGFR bound to gefitinib (Fig. 3, red dots). Representative structures of local energy minima (Fig. 3: 1, 3, and 4) and a saddle-point dividing active and inactive states (Fig. 3: 2) are shown; the major structural rearrangement between the active and the inactive states involves the αC helix and the activation loop region, shown earlier in Fig. 2. The result is in agreement with a previous simulation that



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: Conformational free energy landscape and representative snapshots of WT EGFR computed by metadynamics MD simulations. RM SD is in ˚ A and energy unit is kcal/mol. Contour levels are drawn every 2 kcal/mol. Proteins are rendered using VMD software. 45 showed the Src/Cdk-like inactive state is the global free energy minimum of WT EGFR. 25


Page 10 of 44

Page 11 of 44

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 4: Free energy profile of EGFR interacting with gefitinib.

2.3

Relationship of conformational energy difference and intrinsic binding energies of gefitinib to EGFR WT, L858R and T790M mutants

Fig. 4 combines the energetic results of the TI simulations and the metadynamics MD simulations of EGFRs bound to gefitinib. Sutto and Gervasio demonstrated that the T790M point mutation causes constitutive activation. The mutation stabilizes the active conformation by enhancing the hydrophobic contact between N-lobe and αC helix. 25 The relative energies of the inactive states of WT and the mutants are shown to the left. Upon introduction of the T790M mutation, the active state of EGFR is lower in free energy than the Src/Cdk-like inactive state (Fig. 4, green): the active state and an intermediate conformation with disordered αC helix are isoenergetic. Additional 200 ns constant temperature MD simulation of the T790M mutant bound to gefitinib shows a variation in the conformation of αC-helix while the activation loop keeps its conformation similar to that of the active state (Fig. S7). This result shows that gefitinib can interact with the two isoenergetic lowest



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 states found by Sutto and Gervasio. We expect the binding kinetics of gefitinib to the T790M is diffusion limited. This is markedly different from the binding kinetics to the WT, where the conformational change from the energetically favorable inactive conformation to the unfavorable active conformation may play an important role. Sutto and Gervasio also showed that the active state of the L858R mutant is less favored by 2 kcal/mol relative to semiclosed conformation (Fig. 4, blue). The semiclosed conformation has a markedly distinct conformation of the activation loop compared to that of the active state that can bind to gefitinib. 25 In Fig. 4, the energies of the active conformational states of WT, the L858R and the T790M mutants are placed at 0 for convenience. The levels on the far right converges to the free energy of binding of gefitinib to the active state of each protein. We compute the binding free energy to the active conformation, by combining free energy differences presented here and experimental binding energy (∆Gbind,exp ) in Table 1. We predict the binding free energies to the active state of WT (L858R) and the T790M mutant are −15 and −12 kcal/mol, respectively. Table 3: Binding free energy of gefitinib. EGFR ∆GI→A WT 5 0c T790M L858R 2c

∆Gbind,calcd −10 −12 −13

Kd a ∆Gbind,exp 35.3 −10.2 4.6 −11.4 2.4 −11.8

b

∆∆Gbind,exp M U T ←W T 0 −1.2 −1.6

a

Kd is dissociation constant of gefitinib reported by Yun et al. 17 ∆Gbind = −RT log(Kd ). c Based on metadynamics MD simulation results of Sutto and Gervasio. 25 b

Table 3 compares the computed and the experimental binding affinities. WT is inactive in the apo state, and the overall binding energy of gefitinib is −10 kcal/mol. The T790M mutant is constitutively activated and exists as the active form in apo. The computed binding energy of gefitinib to the T790M mutant is −12 kcal/mol. These quantities can be compared to experimental binding energies. We predict gefitinib binds more favorably to T790M than 12 ACS Paragon Plus Environment

Page 12 of 44

Page 13 of 44

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


to WT by about 2 kcal/mol. Kd values of gefitinib to WT and T790M corresponds to −10.2 and −11.4 kcal/mol at 300 K, respectively. The difference, −1.2 kcal/mol, is close to 2 kcal/mol computed here. TI simulations mutating L858 to R858 using an active state EGFR crystal structure cocrystalized with gefitinib (PDB code: 2ITZ) predicts ∆∆GActive L858R = 0.2 ± 0.9 kcal/mol (Table 2), and the binding free energy is predicted as −13 kcal/mol (Table 3). The reduction in Kd is consistent with the fact that gefitinib sensitivity increases with L858R mutant. 40,43 Thus constitutive activation is the reason that the L858R mutant has lower dissociation constant (Kd ) than WT. The results presented here are relevant to the proposition that the decreased binding affinity is the cause of acquired drug resistance caused by T790M mutation. 8 We find instead that while the T790M single mutation increases the binding affinity of gefitinib to the active form of the protein by 3 kcal/mol, the net binding energy is −12 kcal/mol relative to that of the resting (active) state of the T790M mutation, because of the change in resting state upon mutation. The binding of gefitinib to the T790M mutant is actually 2 kcal/mol more favorable than binding to the WT. 17

2.4

The conformational change of the αC-helix and the binding energy difference explain the binding energy profile of AEE788

Unlike the increased binding effect on gefitinib, the T790M mutation reduces the binding affinity of AEE788 (Table 1). Here we demonstrate that the binding of AEE788 causes a conformational change in the αC-helix, that results in an additional 2 kcal/mol unfavorable energy to the binding energy of AEE788. The energy penalty and the steric effect explain the reduced binding affinity of AEE788 to the T790M mutant relative to the WT. Fig. 5 (a) compares the cocrystal structures of the WT and the T790M mutant EGFR



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

Page 14 of 44

Figure 5: Comparison of WT and T790M EGFR bound to AEE788. (a) Cocrystal structures of the WT (PDB code: 2J6M) and T790M mutant (PDB code: 2JIU) bound to AEE are overlayed. Spheres represent Cα atoms of K745 and I759. (b) Cα-Cα distances between K745 and I759 measured from constant temperature MD simulations. (c) Free energy profile of the Cα-Cα distance measured from metadynamics MD simulations. bound to AEE788. The backbone RMS deviation between the two structures is only 0.87 ˚ A. Also the activation loop (yellow) of the two structures overlay well with each other. However, the αC-helix of the T790M bound to AEE788 tilts away from the inhibitor binding pocket relative to the WT bound to AEE788. We quantify the relative orientation of the αC-helix by measuring the distance between the Cα carbons of K745 and I759: ~ ~ D(K745, I759) = Cα − Cα K745 I759

(3)

~ denotes the 3-dimensional coordinates of the Cα atoms. The distances from the , where Cα WT and the T790M mutant crystal structures are 9 ˚ A and 12 ˚ A, respectively. Fig. 5 (b) summarizes the constant temperature MD simulations of the WT and the


Page 15 of 44

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


T790M bound to AEE788. The orientation of αC-helices observed from the WT cocrystal structure remains stable up to 250 ns (red). When averaged over the last 100 ns, D(K745, I759) of the WT is 9 ˚ A. On the other hand, D(K745, I759) of the T790M mutant bound to AEE788 was extended to 13.6 ˚ A (green), when averaged over the last 100 ns. The conformational change of the αC-helix occurs spontaneously upon binding to the T790M mutant. Additional MD simulations of the T790M bound to AEE788 were carried out using the WT cocrystal structure (PDB code: 2J6M) as the starting conformation (blue). The Cα-Cα distance (D(K745, I759)) increased from 9 ˚ A to > 12 ˚ A in less than 5 ns, that shows the αC-helix spontaneously rotates away from the N-lobe of the protein. These findings confirmed that the conformational change in the αC-helix is necessary for the binding of AEE788 to the T790M mutant EGFR. We performed metadynamics MD simulations to quantify the binding energy change of AEE788 caused by the difference in the conformation of the αC-helix (Fig. 5 (c)). Using the cocrystal structure of WT EGFR bound to AEE788, the Cα atoms distance between K745 and I759 was varied from 7 to 14 ˚ A for 300 ns. The Gibbs free energy profile of D(K745, I759) shows a global minimum at 9 ˚ A that agrees well with the constant temperature MD simulations of WT EGFR bound to AEE788. The free energy profile has additional minimum at 12.7 ˚ A, which coincides with the distance sampled from the constant temperature MD simulations of the T790M mutant bound to AEE788. The difference in free energy between the two minima is 2 kcal/mol, which explains the energy penalty associated with the tilting of the αC-helix. Table. 4 summarized the binding energy profile of AEE788 to WT, T790M, and L858R mutant EGFRs. Yun et al. showed that the dissociation constant (Kd ) of AEE788 to the WT is 5.3 nM, which is ∆Gbind,exp = −11.3 kcal/mol. We predict the binding energy of AEE788 to the WT and the T790M mutant is −11 kcal/mol. The T790M mutation causes 5 kcal/mol energy penalty: 2 kcal/mol due to the conformational change of the αC-helix and 3 kcal/mol due to the binding energy difference to the T790M mutant (Table 2). The



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 16 of 44

Table 4: Binding free energy of AEE788. EGFR ∆GI→A WT 5 0c T790M L858R 2c

∆Gbind,calcd −11 −11 −14

Kd a ∆Gbind,exp 5.3 −11.3 27.3 −10.3 1.1 −12.2

b

∆∆Gbind,exp M U T ←W T 0 +1.0 −0.9

a

Kd [nM] is dissociation constant of AEE788 reported by Yun et al. 17 ∆Gbind = −RT log(Kd ) in kcal/mol. c Based on metadynamics MD simulation results of Sutto and Gervasio. 25 b

gross energy penalty compensates the effect of the constitutive activation caused by the T790M mutation that lowers the binding energy by 5 kcal/mol (Table 3). The experimental dissociation constant of AEE788 to the T790M mutant is 27.3 nM, which corresponds to 1 kcal/mol higher binding energy than to the WT (Table 1). Similarly to gefitinib, we predict the influence of L858R mutation on the binding energy to the active state is negligible. We computed ∆∆GActive L858R = 0.0 ± 1.1 kcal/mol based on TI simulations using the L858R mutant cocrystalized with AEE788 (PDB code: 2JIU). The binding energy of AEE788 to the L858R mutant is lower by 3 kcal/mol due to the constitutive activation. The prediction is in qualitative agreement with the observed Kd profile of AEE788, which shows the binding affinity to the L858R mutant is 1 kcal/mol lower than that of the WT.

2.5

Gatekeeper residue conformation and interaction with inhibitor influences inhibitor binding

Here we detail simulations that show why the active conformation of T790M binds gefitinib 3 kcal/mol less tightly than the WT. Constant temperature MD simulations were performed for EGFR with a bound gefitinib (Fig. 6), which shows the gatekeeper residues adopt different χ1 dihedral angle in the presence of gefitinib (Fig. 7). The χ1 dihedral angle of the threonine gatekeeper prefers gauche+ in apo, whereas trans is favored upon binding of gefitinib. The χ1 angle of the T790M mutant shows different preference. In apo, gauche+/gauche- angles are sampled, whereas all three 16 ACS Paragon Plus Environment

Page 17 of 44

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 6: χ1 dihedral dynamics from constant temperature MD simulations of WT EGFR and T790M mutant: (a) WT EGFR in apo; (b) WT EGFR bound to gefitinib; (c) T790M mutant in apo; (d) T790M mutant bound to gefitinib. conformations are observed upon binding of gefitinib. The finding is in agreement with available cocrystal structures of WT EGFR bound to gefitinib. (Fig. 8). For WT apo structure in the active state, (Fig. 8 a) the χ1 rotamer is gauche+. In the presence of gefitinib (Fig. 8 b), the χ1 dihedral angle of the threonine gatekeeper rotate to trans conformations. For the methionine gatekeeper (M790), the χ1 rotamers in the crystal structures are in gauche- conformations in the apo state as well as in the complex with gefitinib (Fig. 8 c and d). The comparison of MD simulation and the cocrystal structures of the T790M mutant indicates that the preferred χ1 dihedral angle (gauche-) is less stable upon binding of gefitinib. Table 5 summarizes the relative free energies of the gatekeeper sidechain conformations obtained from metadynamics simulations of apo and gefitinib bound complexes. For the case of WT, the preferred χ1 angle of T790 is gauche+, as was found in the with apo state



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

Figure 7: Newman projection of the χ1 dihedral angles of threonine and methionine. crystal structure of EGFR. The trans conformation is 0.7 kcal/mol higher than gauche+, and the gauche- conformation is 2.9 kcal/mol higher than gauche+. Two of the preferred χ1 sidechain dihedral angles were preserved in the presence of gefitinib. The most favorable χ1 is trans as in the gefitinib-bound cocrystal structure of EGFR (PDB code: 2ITO). Gauche+, the lowest in apo form of the WT EGFR, is 0.9 kcal/mol above the trans conformation. The methionine gatekeeper showed a different pattern from the threonine (Table 5). In apo, the trans conformation is 3.9 kcal/mol higher than the gauche+ conformations. However trans conformation becomes isoenergetic to gauche- conformations in the presence of gefitinib, and 0.9–3.9 kcal/mol destabilization results from the sidechain conformational change. This indicates that the preferred gauche+ conformation is forced to the trans to bind gefitinib. Similarly, when AEE788 is bound, which has an aromatic group interacting with the gatekeeper residues, the trans conformation of methionine gatekeeper is isoenergetic with gauche+/gauche-. 18 ACS Paragon Plus Environment

Page 18 of 44

Page 19 of 44

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 8: Conformational changes of the gatekeeper residues upon binding of gefitinib: (a) threonine gatekeeper in apo (PDB Code: 2EB2), (b) threonine gatekeeper in the presence of gefitinib (PDB Code: 2ITO), (c) methionine gatekeeper in apo (PDB Code: 3UG1), and (d) methionine gatekeeper in the presence of gefitinib (PDB Code: 3UG2).

By contrast, binding of ATP does not alter the sidechain conformational preference as in Table 5. With or without ATP, the gauche+ rotamer of both the threonine and the methionine gatekeeper is favored. The computed ∆∆GActive T 790M is close to zero (Table. 2), which shows weaker interaction between the gatekeeper residues and ATP. We also compared the sidechain rotamer preference of the gatekeeper residues with WZ4002 and CNX2006. 46 These compounds were originally designed to target the T790M oncogenic mutation, showing increased potency to EGFRs of both threonine and methionine mutants. 38 As shown in Fig. S10, these inhibitors induce only minor changes in the energy profiles of χ1 dihedral angle with either gatekeepers.



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

Table 5: Energies of different conformations of the gatekeeper residues.

WT apo a T790M apo a WT+Gefitinib b T790M+Gefitinib b WT+AEE788 c T790M+AEE788 c WT+ATP d T790M+ATP d

χ1 rotamer [kcal/mol] Trans Gauche- Gauche+ 0.7 2.9 0.0 3.9 0.9 0.0 0.0 3.8 1.4 0.1 0.0 0.9 0.0 3.1 1.5 0.0 0.5 0.8 0.6 1.0 0.0 2.5 0.7 0.0

a

A gefitinib bound cocrystal structure without the inhibitor was simulated (PDB code: 2ITO). b A gefitinib bound cocrystal structure was simulated (PDB code: 2ITO). c An AEE788 bound cocrystal structure was simulated (PDB code: 2J6M). d An ATP bound cocrystal structure was simulated (PDB code: 2EB3).

2.6

The haloaniline substituent of gefitinib causes unfavorable conformational change of the methionine gatekeeper

The threonine-selective inhibitors (gefitinib and AEE788) induce conformational changes of the sidechains of the gatekeeper residues. Here we analyze the origins of these changes. Upon binding to EGFR, the bulky aromatic group, 3-chloro-4-fluoroaniline of gefitinib, or the phenylethylnamine of AEE788, is placed in close proximity to the gatekeeper. In order to quantify the interaction between the bulky aromatic groups and the gatekeeper residues, we calculated the energetic contribution of the haloaniline substituent of gefitinib (∆∆GActive Haloaniline ) to three possible conformational states of each gatekeeper residue (Figs. 9 and S9). A series of TI simulations were performed in which the haloaniline group was grown onto the truncated gefitinib by turning on the vdW interactions and then the electrostatic interactions of the moiety. The computations were repeated for three different conformational states of the gatekeeper residues, by constraining χ1 angles of the gatekeeper residues to one 20 ACS Paragon Plus Environment

Page 20 of 44

Page 21 of 44

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 9: TI simulation of 3-chloro-4-fluoroaniline substitution. of three preferred χ1 angles (Fig. 7) using a harmonic potential

U (χ1 ) =

k (χ1 − χ∗1 )2 , 2

(4)

χ∗1 was −60◦ , 60◦ , and 180◦ for gauche-, gauche+, and trans conformations respectively, and k was 80 kcal/mol/deg2 . An EGFR structure bound to gefitinib (PDB code: 2ITO) was chosen for the TI simulations. The free energy difference was decomposed into two contributions ∆∆GActive Haloaniline = ∆∆GElec + ∆∆GV dW ,

(5)

∆∆GElec is the electrostatic and ∆∆GV dW is the van der Waals contribution to the overall free energy change. Table 6 summarizes the TI simulations. The haloaniline group contributes favorably to the binding to the threonine gatekeeper by ∆∆GActive Haloaniline = −10 to −7 kcal/mol. The most favorable interaction (∆∆GActive Haloaniline = −9.6 kcal/mol) occurs with the trans rotamer of the threonine. This finding explains why the trans rotamer is the preferred rotamer for the binding of gefitinib. In the apo state, the trans rotamer is only 0.7 kcal/mol higher in 21 ACS Paragon Plus Environment


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

Page 22 of 44

Table 6: Components of interaction energies from TI for creation of of 3-chloro-4fluoroanilinyl group of gefitinib.

∆∆GHaloaniline ∆∆GElec ∆∆GV dW

Trans -9.6 -0.4 -9.2

T790 Gauche- Gauche+ -8.2 -7.2 -0.8 -0.6 -7.4 -6.6

Trans -8.6 -1.0 -7.6

T790M Gauche- Gauche+ -5.6 -4.3 -0.8 -1.1 -4.8 -3.2

energy than gauche+, but the energy penalty is compensated by more favorable interaction with the trans rotamer. For the methionine gatekeeper, the most favorable interaction also occurs with the trans rotamer. The magnitude of the interaction (−8.6 kcal/mol) was comparable to that of the threonine gatekeeper, but only at the expense of the 3.9 kcal/mol destabilization penalty to achieve the conformation (Table 5). Gauche- and gauche+ rotamers forced the methanethiyl group of the methionine sidechain to rotate toward the haloaniline of gefitinib, which causes unfavorable steric effects. Thus none of the three possible sidechain conformations confer favorable interaction with the haloaniline group of gefitinib, which explains unfavorable steric effect (∆∆Gbind,Active = 2–4 kcal/mol) caused by the gatekeeper mutation (Table 3). T 790M The decomposition analysis shows the van der Waals interaction (∆∆GV dW ) contributes the most, whereas the electrostatic (∆∆GElec ) was almost identical regardless of the gatekeeper composition or the χ1 rotamer conformation.

2.7

The origin of the high energy of the trans rotamer of the methionine gatekeeper

In Table 7, we compared the conformational preferences of the methionine gatekeeper computed from metadynamics simulations to that of known protein structures using sidechain conformational preference (rotamer library) compiled by Shapovalov and Dunbrack. 47 The sidechain conformational preference is defined as a probability of finding a quartet of sidechain


Page 23 of 44

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


Table 7: Rotamer energies from the sidechains of the gatekeeper residues from metadynamics simulations of EGFR and from the Dunbrack library of rotamer preference Trans Metadynamics 0.7 Dunbrack 1.2

T790 Gauche- Gauche+ 2.9 0.0 2.4 0.0

Trans 3.9 0.9

M790 Gauche- Gauche+ 0.9 0.0 0.0 0.6

dihedral angles (χ1 , χ2 , χ3 , and χ4 ) for a given set of backbone torsional angles (φ, ψ)

PDunbrack (χ1 , χ2 , χ3 , χ4 |φ, ψ).

(6)

The free energy of the χ1 dihedral angle was computed at 300 K from the Dunbrack’s results

P (χ1, χ2, χ3 χ4) =

X

PDunbrack (χ1 , χ2 , χ3 , χ4 |φ, ψ)

(7)

φ,ψ

! E(χ1 ) = −RT log

X

P (χ1 , χ2 , χ3 , χ4 )

(8)

χ2 ,χ3 ,χ4

, where P (χ1, χ2, χ3, χ4) is probability of finding a certain rotamer quartet. The summation in Eq. 7 was computed over the pairs of φ and ψ angles sampled from 60 ns constanttemperature MD simulations. The conformational energy of the trans rotamer of methionine is 3.0 kcal/mol higher in T790M of EGFR than in available high-resolution crystal structures. Additional MD simulations of T790M were performed with an added restraining potential to keep the methionine gatekeeper in the the trans χ1 dihedral angle using Eq. 4 with the spring constant, k, was 80 kcal/mol/deg2 . The analysis of the MD trajectory revealed the presence of steric clash between the Cγ-H atom of the methionine and backbone N H a neighboring glutamine (Q791): the distance between a pair of atoms ( Cγ of M790: and N of Q791 ) stays within 2.0 ˚ A over 30 ns, and in 1.7 ˚ A in the snapshot shown in Fig. 10.



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 10: The trans rotamer of the methionine gatekeeper experiences steric clash.

2.8

T790M mutant stabilizes hydrophobic contact with the active state

We surveyed how the T790M mutation causes the conformational equilibrium in favor of the active state inapo. Fig. 11 (a) describes the αC-helix, the hydrophobic core, the activation loop, and the phosphate binding loop. 48 Both the αC-helix and the activation loop are explained in introduction. The hydrophobic core mediates contacts between N-lobe and αChelix. Upon binding of ATP, the phosphate binding loop makes contact with the phosphate group. Fig. S4 compares the conformational fluctuation of the selected regions of the WT and the T790M mutant. The hydrophobic core of the T790M mutant deviates less from the initial structure of the MD simulations than the WT: RMS deviation of the T790M is 2.7 ˚ A, whereas RMSD of the WT is 3.4 ˚ A, when averaged over 200 ns. On the other hand, the average RMS deviations of the αC-helix, the activation loop, and the p-loop of both the WT and the T790M mutant were comparable. The finding suggests that the T790M mutant stabilizes the hydrophobic core. Also when MD simulation snapshots are overlayed, the N-terminal fragment of αC-helix of the WT is disordered relative to that of the T790M


Page 24 of 44

Page 25 of 44

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 11: (a) Structural elements of EGFR. We adopted amino acid sequences defined by Shih et al.: αC-helix, P757-A767; hydrophobic core, L747, E749, T751, K745, A763, M766, A767, V769, L777, I780, L788, G857, L858, L861, L862; activation loop, I155-I180; phosphate binding loop (G-loop), G719-G724. 48 (b) Overlay of αC-helix of WT (magenta) and T790M (cyan) EGFR sampled from 250 ns MD simulation trajectories. mutant, as shown in Fig. 11 (b). These findings are in line with a previous metadynamics MD simulation result, where increased stability of the hydrophobic patch between N-lobe and αC-helix explains the increased stability of the active state of T790M. 25

3

Conclusions

Using thermodynamic integration (TI) MD simulations, we show that the gatekeeper mutation (T790M) causes a 2–3 kcal/mol energy penalty for binding of certain EGFR inhibitors to the active conformation of EGFR. The gatekeeper mutation causes the constitutive activation, shifting the conformational equilibrium toward the active state of EGFR. The steric



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

effect and the constitutive activation together explain observed binding affinity profiles of gefitinib to WT, and T790M and L858R mutants. MD simulations show that the T790M mutation causes a conformational change of the αC-helix upon binding of AEE788. The steric effect and the energy penalty of deforming the αC-helix explain unfavorable binding of AEE788 to the T790M mutant compared to the WT. We also demonstrate that upon binding of an EGFR inhibitor with an aromatic group that makes contact with the gatekeeper residue, the sidechain conformation of the WT and the mutated gatekeeper changes. The conformational changes in the gatekeeper residues are analyzed using TI and metadynamics MD simulations. The χ1 dihedral angle of the threonine gatekeeper can change to an energetically favored angle. However the χ1 dihedral angle of the methionine gatekeeper shifts to an energetically unfavorable value. The repulsion between the sidechain of the methionine gatekeeper and the aromatic substituents of EGFR inhibitors cause the energetically unfavorable conformational change. Finally we show that the T790M stabilizes the hydrophobic contacts between the substrate binding site and the αC-helix of the active state. The enhanced stability of the hydrophobic contacts explains the constitutive activation of the active state by the gatekeeper mutation. Chronic kinase inhibition in oncology treatment engenders the emergence of tumor cells expressing kinase mutants that circumvent therapeutic inhibition. The influence of mutations on kinase function, and the therapeutic effectiveness of kinase inhibitors can be traced to a number of important factors. By dissecting the physical consequences of mutation on the observed Kd of inhibitors in the wild-type and the mutant EGFR isoforms, we hope to have identified novel design strategies to thwart mutational resistance. Specifically, we have observed significant alterations in conformational preferences for the target protein, and downstream consequences to substrate and inhibitor affinity upon gatekeeper mutation. Upon examination of the energetics leading to the influence of mutations on inhibitor potency, we have shown that conformational equilibrium shift toward the active state, steric effects on the inhibitor-protein binding, and gatekeeper sidechain conformational changes


Page 26 of 44

Page 27 of 44

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


induced by the inhibitor binding all contribute to the binding energy change in a significant way. In traditional inhibitor design, there is often an emphasis on one or two of these factors, but where mutation impacts drug design, the subtle interplay of all of these factors is important.

4

Methods

4.1

Preparation of MD simulations

Crystal structures of active-state EGFR kinase were retrieved from PDB database for the MD simulations (Table S1). The ff99SB 49 force fields was used for the EGFR structures. Gatekeeper mutations were introduced using tleap program of Amber14 50 suite. The generalized amber force field (GAFF) 51 was used for EGFR inhibitors. Atomic partial charges of the inhibitors were assigned based on the RESP model. 52 Starting from the crystallographic coordinates of each inhibitor-protein complexes, optimization of the inhibitor geometry was performed with density functional theory at the B3LYP/6-31G*; single point energies and charges were computed with HF/6-31G* theory. We used Gaussian09 53 for the quantum mechanical (QM) computations and antechamber 54 module in Amber14 to process atomic partial charges from the QM computations. For the thermodynamic integration (TI) simulations, apo and protein-ligand complexes (Table. S1) were solvated in octahedral boxes filled with TIP3P water molecules. We ensured 11 ˚ A margin from the protein surface to the solvation boundary. Counter ions (Cl- or Na+) were added to neutralize the solvated molecules. An active state (PDB code: 2ITO) and an inactive state (PDB code: 2GS7) crystallographic structures of EGFR were used as the initial models for the constant temperature MD simulations and the metadynamics MD simulations. For apo simulations, bound ligands in the crystal structures were removed. Residues from P699 to G982 were used to make the amino acid sequences of the two states are identical. Missing residues in the crystallo27 ACS Paragon Plus Environment


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

Page 28 of 44

graphic structures were built using ModLoop web server. 55 For the studies of AEE788, two of cocrystal structures of WT and T790M mutant were used (PDB codes: 2J6M and 2JIU). Of note, fractions of the activation loop residues are missing from the crystallographic structures: G866 to G876 of the active state and H870 to K875 of the inactive state. Long time MD simulations done by Shan et al. showed that the activation loop is disordered. 24 However the MD trajectories initiated from a specific conformation may not sample entire conformational ensemble of the activation loop. We compared the conformational ensemble of the activation loop sampled from additional two MD trajectories. The missing residues of the activation loop are remodeled by using two snapshots from constant temperature MD simulations as the templates for ModLoop web server (Fig. S2). 200 ns MD simulations were carried out starting from the newly generated models. As shown in Fig. S3, the repeated MD simulations sampled comparable region in the RMSD plane, showing that the initial loop model does not cause a bias in MD simulations.

4.2

Thermodynamic integration (TI) simulations

We used the thermodynamic integration (TI) 30,56 method to compute the binding free energy difference to the active-state of EGFR upon gatekeeper mutation (∆∆GActive T 790M ) or different sidechain conformation of the gatekeeper residue (∆∆GActive Haloaniline ). The TI module in Amber14 57 was used. As shown in Fig. S1, we simulate alchemical transformation from an initial state (i.e. T790) to a final state (i.e. M790), using hybrid Hamiltonians; 30

H(λ) = (1 − λ)HT 790 + λHM 790 (λ ∈ {λ|0 . . . 1})

(9)

, where HT 790 and HM 790 stands for the solvated Hamiltonians of the WT and the gatekeeper mutant proteins, and λ was varied from 0 to 1, which defines hypothetical intermediate states of the alchemical transformation. We simulated two alchemical transformations for


Page 29 of 44

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


each problem: (1) without inhibitor (∆G1 ), and with an inhibitor (∆G2 ) for the studies of ∆∆GActive T 790M (Sec. 2.1); (2) without or with the 3-chloro-4-fluoro-aniline moiety both within the ligand binding pocket and within the solvent (Sec. 2.6). The relative binding free energy is defined: ∆∆G = ∆G2 − ∆G1

(10)

, where ∆G1 and ∆G2 are initial and final state of the alchemical transformation, respectively. For each alchemical transformation, electrostatic and vdW contributions were decoupled into three individual steps:

HR Step 1 ∆GTElec : Turning-off of electrostatic interactions of the initial state ET Step 2 ∆GTV HR→M : VdW transformation of the initial and the final state dW ET Step 3 ∆GM Elec : Turning-on of electrostatic interactions of the final state

The decoupling of electrostatic and vdW interaction was used in order to improve statistical convergence of TI simulations. 58 For Steps 1 and 3, λ was sampled in 0.05 increment from 0 to 1. The vdW transformation suffers more statistical uncertainty, 56 so that λ was sampled in 0.04 increment and softcore vdW potential 59 was used for Step 2. The decoupling scheme leads to the free energy decomposition analysis:

∆∆GElec = ∆∆GV dW =

HR ET T HR M ET ∆GTElec + ∆GM Elec ∆G2 − ∆GElec + ∆GElec ∆G1 ET T HR→M ET ∆GTV HR→M − ∆G dW V dW ∆G2 ∆G1

(11) (12)

, where ∆∆GElec is the difference of electrostatic interactions between the final and the initial state in the presence of vdW interaction and ∆∆GV dW is the difference of the vdW interaction in the absence of electrostatic interaction.



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 30 of 44

The free energy difference of the end states was computed; Z

∂H ∆G = dλ ∂λ λ X ∂H ' 4λ, ∂λ λ

(13) (14)

, where 4λ is the increment between successive λ values: 0.05 for Steps 1 and 3, and 0.04 for Step 2. The numerical quadrature (Eq. 14) was performed with the trapezoidal rule. 60 A more sophisticated quadrature (Simpson’s rule) was tested but gave identical results. For the analysis of statistical uncertainty, we followed the running average procedure. 61 For each λ point, a running average of h∂H/∂λiλ on every 50 ps interval was measured. The variance of the running average was understood as the uncertainty of each λ point: !2 N X 1 ∂H σλ2 = − N i=1 ∂λ λ,i

N 1 X ∂H N i=1 ∂λ λ,i

!2 (15)

, where i = 1, . . . , N is the index of each 50 ps interval. The variance (σ 2 ) of the computed free energy difference was computed;

σ2 =

X

σλ2 (4λ)2

(16)

, as each λ window is statistically independent. During the course of TI simulations, we monitored ∆∆G computed at time t is within ±σ/2 of the ∆∆G computed at t-1 ns (Fig. S11). Steps 1 and 3 used to converge withing 3 ns. However step2 converges slower than the others, that sometimes the simulation time increased to 6 ns. Before initiating TI simulations, a solvated system was energy minimized in 4000 steps. All bonds connecting hydrogens and heavy atoms were restrained using the SHAKE algorithm. Each system was equilibrated for 1 ns at 300 K while restraining Cα atoms to their 2

crystallographic coordinates using 1.0 kcal/mol/˚ A spring constant. Restraints on backbone atoms reduces conformational flexibility of the protein and thus lowers statistical uncer30 ACS Paragon Plus Environment

Page 31 of 44

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


tainty of the TI calculations. 36 The Langevin thermostat with the 5.0 ps collision frequency was used to control the simulation temperature at 300 K. A subsequent 1 ns equilibrations was followed with Berendsen barostat 61 to control the pressure at 1 atm, with weaker re2

straints on Cα atoms: spring constant was 0.2 kcal/mol/˚ A . After the equilibration period, h∂H/∂λiλ was computed every 1 ps, which lasted 1 to 5 ns. Amber14 was used for the TI simulations.

4.3

Metadynamics simulations of conformational equilibrium of EGFR

Metadynamics MD simulations 29 were carried out to quantify the free energy differences of the conformational states of EGFR. NAMD v2.9 62 was used for the metadynamics simulations with ff99SB force fields. 49 Well-tempered metadynamics 63 in conjunction with multiplewalker algorithm 64 were used to accelerate convergence of the free energy landscape computed from the metadynamics simulations. The Langevin thermostat controlled the temperature at 300 K and Langevin barostat kept the pressure at 1 atm. The free energy difference of the active and the Src/Cdk-like inactive states is computed:

∆GI→A = GActive − GInactive

(17)

, where GActive and GInactive are conformational free energy of the active and the Src/Cdk-like inactive state of EGFR, respectively (Sec. 2.2). We simulated 18 replicas of EGFR, where half them started from the active state and the rest began from the inactive state. We chose three collective variables to describe the conformational equilibrium of EGFR:

RM SD1 : RMS deviation to the active state RM SD2 : RMS deviation to the inactive state dP C : Dihedral-angle principle component. 31 ACS Paragon Plus Environment


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

The range of the collective variables were constrained, to facilitate the conformational transition:

−3.2 ˚ A ≤ (RM SD1 − RM SD2) ≤ 3.2 ˚ A 4.5 ˚ A ≤ (RM SD1 + RM SD2) ≤ 7.0 ˚ A −3.0 ≤

dP C

≤ 3.0

Amino acids within 5.0 ˚ A from the ligand binding pocket or the region forming α-helix and β-sheet secondary structure were controlled by the restraints provided by metadynamics algorithm (Table S3). The dihedral angle principle component (dP C) 65,66 was introduced to accelerated conformational transition of the activation-loop region. In the active state, the beginning of the activation-loop (G858-L863) is unstructured, whereas in the inactive state the same region forms α-helix (Fig. 2). Initial metadynamics simulations using RM SD based collective variables failed to trigger the conformation transition of the activation-loop. We recognize dP C complements RM SD1 and RM SD2 to accelerate further the conformational transition between the active and the inactive states of EGFR. We used grcarma software 67 for the dihedral principle component analysis of the active and the inactive state of EGFR. After 400 steps of energy minimization, each solvated system was equilibrated in 12 ns before turning on metadynamics simulation which lasted up to 150 ns per replica. Metadynamics MD simulations were used to quantify the free energy difference of the distinct αC-helix conformations in the presence of AEE788 (Sec. 2.4). The Cα-Cα distance of K745 and I759 was selected as the reaction coordinate describing the αC-helix conformation: ~ ~ D(K745, I759) = Cα − Cα K745 I759 . Before initiating metadynamics MD simulations, a crystallographic structure of WT EGFR cocrystalized with AEE788 (PDB code: 2J6M) was heated and equilibrated in 20 ns. We simulated 6 replicas of the EGFR-AEE788 complex for 180 ns in total. 32 ACS Paragon Plus Environment

Page 32 of 44

Page 33 of 44

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


In order to study the influences of the sidechain conformation adaptation upon binding of inhibitors (Sec. 2.5), the χ1 dihedral angles of gatekeeper residues were chosen as the collective variable of the metadynamics simulations. Four consecutive atoms of threonine and methionine define χ1 dihedral angle:

χ1 of threonine :

N –Cα–Cβ–Oγ

χ1 of methionine : N –Cα–Cβ–Cγ

Both apo and inhibitor bound complex were simulated. Each solvated system was energy minimized in 400 steps, followed by 1 ns equilibration with restraints on Cα atoms to their crystallographic coordinates. Additional 2 ns equilibration without restraints were then followed. After completion of the equilibration period, metadynamics MD simulations were carried out for 100 ns per each replica.

4.4

Computation of the binding energy of gefitinib and AEE788

Binding free energy of an inhibitor (∆Gbind ) is decomposed into three contributions:

∆Gbind = ∆GI→A + ∆Gbind,Active + ∆∆Gbind,Active x x x WT

(18)

, where x is WT or the gatekeeper mutation (T790M) of EGFR, ∆GI→A is the free energy of transforming the inactive state to the active state, ∆Gbind,Active is the binding free energy to WT the active state without gatekeeper mutation, and ∆∆Gbind,Active is the free energy penalty x caused by the mutation. By subtracting the binding free energy to WT (∆Gbind W T ) from the binding free energy to the gatekeeper mutant (∆Gbind x ), we have the relative binding free energy of the gatekeeper mutant:

I→A ∆∆Gbind + ∆∆Gbind,Active W T →x = ∆∆Gx x


(19)


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

I→A , where ∆∆GI→A − ∆GI→A W T →x = ∆Gx W T is the difference of the free energy to transform the

active state to the inactive state.

5

Acknowledgments

J.P. and K.N.H. are grateful for financial support from Celgene, Bedford, MA. Computational resources, Hoffman2 cluster from UCLA Institute for Digital Research and Education and Keeneland supercomputer through XSEDE program (TG-CHE040013N), are gratefully acknowledged.

6

Supporting Information

A complete list of the PDB coordinates used for the manuscript, the duration and types of MD simulations, and the selected amino acids for the RMSD computations are summarized in the supporting tables. Supporting figures show details of the TI and the metadynamics MD simulations, structural modeling of the missing residues in the PDB coordinates, and constant temperature MD simulations of the T790M mutant. This information is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Herbst, R. S. Review of epidermal growth factor receptor biology. Int. J. Radiat. Oncol. Biol. Phys. 2004, 59, 21–26. (2) Jura, N.; Endres, N. F.; Engel, K.; Deindl, S.; Das, R.; Lamers, M. H.; Wemmer, D. E.; Zhang, X.; Kuriyan, J. Mechanism for activation of the EGF receptor catalytic domain by the juxtamembrane segment. Cell 2009, 137, 1293–1307. (3) Cheng, L.; Alexander, R. E.; Maclennan, G. T.; Cummings, O. W.; Montironi, R.;


Page 34 of 44

Page 35 of 44

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


Lopez-Beltran, A.; Cramer, H. M.; Davidson, D. D.; Zhang, S. Molecular pathology of lung cancer: key to personalized medicine. Mod. Pathol. 2012, 25, 347–369. (4) Morphy, R. Selectively nonselective kinase inhibition: striking the right balance. J. Med. Chem. 2010, 53, 1413–1437. (5) Pao, W.; Miller, V.; Zakowski, M.; Doherty, J.; Politi, K.; Sarkaria, I.; Singh, B.; Heelan, R.; Rusch, V.; Fulton, L.; Mardis, E.; Kupfer, D.; Wilson, R.; Kris, M.; Varmus, H. EGF receptor gene mutations are common in lung cancers from ”never smokers” and are associated with sensitivity of tumors to gefitinib and erlotinib. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 13306–13311. (6) Burris, H. A., 3rd Dual kinase inhibition in the treatment of breast cancer: initial experience with the EGFR/ErbB-2 inhibitor lapatinib. Oncologist 2004, 9 Suppl 3, 10–15. (7) Juchum, M.; G¨ unther, M.; Laufer, S. A. Fighting cancer drug resistance: Opportunities and challenges for mutation-specific EGFR inhibitors. Drug Resist. Updat. 2015, 20, 12–28. (8) Kobayashi, S.; Boggon, T. J.; Dayaram, T.; Jnne, P. A.; Kocher, O.; Meyerson, M.; Johnson, B. E.; Eck, M. J.; Tenen, D. G.; Halmos, B. EGFR mutation and resistance of non-small-cell lung cancer to gefitinib. N. Engl. J. Med. 2005, 352, 786–792. (9) Godin-Heymann, N.; Ulkus, L.; Brannigan, B. W.; McDermott, U.; Lamb, J.; Maheswaran, S.; Settleman, J.; Haber, D. A. The T790M ”gatekeeper” mutation in EGFR mediates resistance to low concentrations of an irreversible EGFR inhibitor. Mol. Cancer Ther. 2008, 7, 874–879. (10) Azam, M.; Seeliger, M. A.; Gray, N. S.; Kuriyan, J.; Daley, G. Q. Activation of tyrosine kinases by mutation of the gatekeeper threonine. Nat. Struct. Mol. Biol. 2008, 15, 1109–1118. 35 ACS Paragon Plus Environment


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

(11) Landi, L.; Cappuzzo, F. Targeted therapies: Front-line therapy in lung cancer with mutations in EGFR. Nat. Rev. Clin. Oncol. 2011, 8, 571–573. (12) Costa, D. B.; Schumer, S. T.; Tenen, D. G.; Kobayashi, S. Differential responses to erlotinib in epidermal growth factor receptor (EGFR)-mutated lung cancers with acquired resistance to gefitinib carrying the L747S or T790M secondary mutations. J. Clin. Oncol. 2008, 26, 1182–4; author reply 1184–6. (13) Brehmer, D.; Greff, Z.; Godl, K.; Blencke, S.; Kurtenbach, A.; Weber, M.; Mller, S.; Klebl, B.; Cotten, M.; Kri, G.; Wissing, J.; Daub, H. Cellular targets of gefitinib. Cancer Res. 2005, 65, 379–382. (14) Köhler, J.; Schuler, M. Afatinib, erlotinib and gefitinib in the first-line therapy of EGFR mutation-positive lung adenocarcinoma: a review. Onkologie 2013, 36, 510–518. (15) Sequist, L. V.; Besse, B.; Lynch, T. J.; Miller, V. A.; Wong, K. K.; Gitlitz, B.; Eaton, K.; Zacharchuk, C.; Freyman, A.; Powell, C.; Ananthakrishnan, R.; Quinn, S.; Soria, J.C. Neratinib, an irreversible pan-ErbB receptor tyrosine kinase inhibitor: results of a phase II trial in patients with advanced non-small-cell lung cancer. J. Clin. Oncol. 2010, 28, 3076–3083. (16) Gajiwala, K. S.; Feng, J.; Ferre, R.; Ryan, K.; Brodsky, O.; Weinrich, S.; Kath, J. C.; Stewart, A. Insights into the aberrant activity of mutant EGFR kinase domain and drug recognition. Structure 2013, 21, 209–219. (17) Yun, C.-H.; Mengwasser, K. E.; Toms, A. V.; Woo, M. S.; Greulich, H.; Wong, K.-K.; Meyerson, M.; Eck, M. J. The T790M mutation in EGFR kinase causes drug resistance by increasing the affinity for ATP. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 2070–2075. (18) Schrödinger, LLC, The PyMOL Molecular Graphics System, Version 1.3r1. 2010.


Page 36 of 44

Page 37 of 44

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


(19) Kumar, A.; Petri, E. T.; Halmos, B.; Boggon, T. J. Structure and clinical relevance of the epidermal growth factor receptor in human cancer. J. Clin. Oncol. 2008, 26, 1742–1751. (20) Tong, M.; Seeliger, M. A. Targeting conformational plasticity of protein kinases. ACS Chem Biol 2015, 10, 190–200. (21) Okabe, T.; Okamoto, I.; Tamura, K.; Terashima, M.; Yoshida, T.; Satoh, T.; Takada, M.; Fukuoka, M.; Nakagawa, K. Differential constitutive activation of the epidermal growth factor receptor in non-small cell lung cancer cells bearing EGFR gene mutation and amplification. Cancer Res 2007, 67, 2046–2053. (22) Chong, C. R.; Jnne, P. A. The quest to overcome resistance to EGFR-targeted therapies in cancer. Nat. Med. 2013, 19, 1389–1400. (23) Zuccotto, F.; Ardini, E.; Casale, E.; Angiolini, M. Through the ”gatekeeper door”: exploiting the active kinase conformation. J. Med. Chem. 2010, 53, 2681–2694. (24) Shan, Y.; Arkhipov, A.; Kim, E. T.; Pan, A. C.; Shaw, D. E. Transitions to catalytically inactive conformations in EGFR kinase. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 7270–7275. (25) Sutto, L.; Gervasio, F. L. Effects of oncogenic mutations on the conformational freeenergy landscape of EGFR kinase. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 10616– 10621. (26) Kim, Y.; Li, Z.; Apetri, M.; Luo, B.; Settleman, J. E.; Anderson, K. S. Temporal resolution of autophosphorylation for normal and oncogenic forms of EGFR and differential effects of gefitinib. Biochemistry 2012, 51, 5212–5222. (27) Carey, K. D.; Garton, A. J.; Romero, M. S.; Kahler, J.; Thomson, S.; Ross, S.; Park, F.; Haley, J. D.; Gibson, N.; Sliwkowski, M. X. Kinetic analysis of epidermal growth factor 37 ACS Paragon Plus Environment


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

receptor somatic mutant proteins shows increased sensitivity to the epidermal growth factor receptor tyrosine kinase inhibitor, erlotinib. Cancer Res. 2006, 66, 8163–8171. (28) Zhang, X.; Pickin, K. A.; Bose, R.; Jura, N.; Cole, P. A.; Kuriyan, J. Inhibition of the EGF receptor by binding of MIG6 to an activating kinase domain interface. Nature 2007, 450, 741–744. (29) Laio, A.; Gervasio, F. L. Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys. 2008, 71, 126601. (30) Kirkwood, J. G. Statistical mechanics of fluid mixtures. J. Chem. Phys. 1935, 3, 300– 313. (31) Chan, W. W.; Wise, S. C.; Kaufman, M. D.; Ahn, Y. M.; Ensinger, C. L.; Haack, T.; Hood, M. M.; Jones, J.; Lord, J. W.; Lu, W. P.; Miller, D.; Patt, W. C.; Smith, B. D.; Petillo, P. A.; Rutkoski, T. J.; Telikepalli, H.; Vogeti, L.; Yao, T.; Chun, L.; Clark, R.; Evangelista, P.; Gavrilescu, L. C.; Lazarides, K.; Zaleskas, V. M.; Stewart, L. J.; Van Etten, R. A.; Flynn, D. L. Conformational control inhibition of the BCR-ABL1 tyrosine kinase, including the gatekeeper T315I mutant, by the switch-control inhibitor DCC2036. Cancer Cell 2011, 19, 556–568. (32) Pedersen, M. W.; Pedersen, N.; Ottesen, L. H.; Poulsen, H. S. Differential response to gefitinib of cells expressing normal EGFR and the mutant EGFRvIII. Br. J. Cancer 2005, 93, 915–923. (33) Li, D.; Ambrogio, L.; Shimamura, T.; Kubo, S.; Takahashi, M.; Chirieac, L. R.; Padera, R. F.; Shapiro, G. I.; Baum, A.; Himmelsbach, F.; Rettig, W. J.; Meyerson, M.; Solca, F.; Greulich, H.; Wong, K.-K. BIBW2992, an irreversible EGFR/HER2 inhibitor highly effective in preclinical lung cancer models. Oncogene 2008, 27, 4702–4711.


Page 38 of 44

Page 39 of 44

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


(34) Traxler, P.; Allegrini, P. R.; Brandt, R.; Brueggen, J.; Cozens, R.; Fabbro, D.; Grosios, K.; Lane, H. A.; McSheehy, P.; Mestan, J.; Meyer, T.; Tang, C.; Wartmann, M.; Wood, J.; Caravatti, G. AEE788: a dual family epidermal growth factor receptor/ErbB2 and vascular endothelial growth factor receptor tyrosine kinase inhibitor with antitumor and antiangiogenic activity. Cancer Res. 2004, 64, 4931–4941. (35) Kancha, R. K.; von Bubnoff, N.; Peschel, C.; Duyster, J. Functional analysis of epidermal growth factor receptor (EGFR) mutations and potential implications for EGFR targeted therapy. Clin. Cancer Res. 2009, 15, 460–467. (36) Zhu, S.; Travis, S. M.; Elcock, A. H. Accurate Calculation of Mutational Effects on the Thermodynamics of Inhibitor Binding to p38α MAP Kinase: A Combined Computational and Experimental Study. J. Chem. Theory Comput. 2013, 9, 3151–3164. (37) Reardon, D. A.; Conrad, C. A.; Cloughesy, T.; Prados, M. D.; Friedman, H. S.; Aldape, K. D.; Mischel, P.; Xia, J.; DiLea, C.; Huang, J.; Mietlowski, W.; Dugan, M.; Chen, W.; Yung, W. K. A. Phase I study of AEE788, a novel multitarget inhibitor of ErbB- and VEGF-receptor-family tyrosine kinases, in recurrent glioblastoma patients. Cancer Chemother. Pharmacol. 2012, 69, 1507–1518. (38) Zhou, W.; Ercan, D.; Chen, L.; Yun, C.-H.; Li, D.; Capelletti, M.; Cortot, A. B.; Chirieac, L.; Iacob, R. E.; Padera, R.; Engen, J. R.; Wong, K.-K.; Eck, M. J.; Gray, N. S.; Jnne, P. A. Novel mutant-selective EGFR kinase inhibitors against EGFR T790M. Nature 2009, 462, 1070–1074. (39) Sakuma, Y.;

Yamazaki, Y.;

Nakamura, Y.;

Yoshihara, M.;

Matsukuma, S.;

Nakayama, H.; Yokose, T.; Kameda, Y.; Koizume, S.; Miyagi, Y. WZ4002, a thirdgeneration EGFR inhibitor, can overcome anoikis resistance in EGFR-mutant lung adenocarcinomas more efficiently than Src inhibitors. Lab. Invest. 2012, 92, 371–383. (40) Gazdar, A. F. Activating and resistance mutations of EGFR in non-small-cell lung 39 ACS Paragon Plus Environment


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

cancer: role in clinical response to EGFR tyrosine kinase inhibitors. Oncogene 2009, 28 Suppl 1, S24–S31. (41) Jorissen, R. N.; Walker, F.; Pouliot, N.; Garrett, T. P. J.; Ward, C. W.; Burgess, A. W. Epidermal growth factor receptor: mechanisms of activation and signalling. Exp. Cell Res. 2003, 284, 31–53. (42) Normanno, N.; De Luca, A.; Bianco, C.; Strizzi, L.; Mancino, M.; Maiello, M. R.; Carotenuto, A.; De Feo, G.; Caponigro, F.; Salomon, D. S. Epidermal growth factor receptor (EGFR) signaling in cancer. Gene 2006, 366, 2–16. (43) Mulloy, R.; Ferrand, A.; Kim, Y.; Sordella, R.; Bell, D. W.; Haber, D. A.; Anderson, K. S.; Settleman, J. Epidermal Growth Factor Receptor Mutants from Human Lung Cancers Exhibit Enhanced Catalytic Activity and Increased Sensitivity to Gefitinib. Cancer Res. 2007, 67, 2325–2330. (44) Gazdar, A. F. Personalized medicine and inhibition of EGFR signaling in lung cancer. N. Engl. J. Med. 2009, 361, 1018–1020. (45) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graph. 1996, 14, 33–38. (46) Ohashi, K.; Suda, K.; Sun, J.; Pan, Y.; Walter, A. O.; Dubrovskiy, A.; Tjin, R.; Mitsudomi, T.; Pao, W. CNX-2006, a novel irreversible epidermal growth factor receptor (EGFR) inhibitor, selectively inhibits EGFR T790M and fails to induce T790Mmediated resistance in vitro. Cancer Research 2013, 73, 2101A. (47) Shapovalov, M. V.; Dunbrack, R. L., Jr A smoothed backbone-dependent rotamer library for proteins derived from adaptive kernel density estimates and regressions. Structure 2011, 19, 844–858.


Page 40 of 44

Page 41 of 44

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


(48) Shih, A. J.; Telesco, S. E.; Choi, S.-H.; Lemmon, M. A.; Radhakrishnan, R. Molecular dynamics analysis of conserved hydrophobic and hydrophilic bond-interaction networks in ErbB family kinases. Biochem J 2011, 436, 241–251. (49) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins 2006, 65, 712–725. (50) Case, D.; Babin, V.; Berryman, J.; Betz, R.; Cai, Q.; Cerutti, D.; T.E. Cheatham, I.; Darden, T.; Duke, R.; Gohlke, H.; Goetz, A.; Gusarov, S.; Homeyer, N.; Janowski, P.; Kaus, J.; Kolossváry, I.; Kovalenko, A.; Lee, T.; LeGrand, S.; Luchko, T.; Luo, R.; Madej, B.; Merz, K.; Paesani, F.; Roe, D.; Roitberg, A.; Sagui, C.; Salomon-Ferrer, R.; Seabra, G.; Simmerling, C.; Smith, W.; Swails, J.; Walker, R.; Wang, J.; Wolf, R.; Wu, X.; Kollman, P. AMBER 14. 2014; University of California, San Francisco. (51) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157–1174. (52) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 1993, 97, 10269–10280. (53) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.;



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

Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian09 Revision D.01. 2009; Gaussian Inc. Wallingford CT. (54) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graph. Model. 2006, 25, 247–260. (55) Fiser, A.; Sali, A. ModLoop: Automated modeling of loops in protein structures. Bioinformatics 2003, 19, 2500–2501. (56) Steinbrecher, T.; Labahn, A. Towards accurate free energy calculations in ligand protein-binding studies. Curr. Med. Chem. 2010, 17, 767–785. (57) Kaus, J. W.; Pierce, L. T.; Walker, R. C.; McCammon, J. A. Improving the Efficiency of Free Energy Calculations in the Amber Molecular Dynamics Package. J. Chem. Theory Comput. 2013, 9, 4131–4139. (58) Wu, K.-W.; Chen, P.-C.; Wang, J.; Sun, Y.-C. Computation of relative binding free energy for an inhibitor and its analogs binding with Erk kinase using thermodynamic integration MD simulation. J. Comput. Aided Mol. Des. 2012, 26, 1159–1169. (59) Zacharias, M.; Straatsma, T. P.; McCammon, J. A. Separation-shifted scaling, a new scaling method for Lennard-Jones interactions in thermodynamic integration. J. Chem. Phys. 1994, 100, 9025–9031. (60) Bruckner, S.; Boresch, S. Efficiency of alchemical free energy simulations. II. Improvements for thermodynamic integration. J. Comput. Chem. 2011, 32, 1320–1333.


Page 42 of 44

Page 43 of 44

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


(61) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: New York, NY, USA, 1989. (62) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kal, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781–1802. (63) Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100, 020603. (64) Raiteri, P.; Laio, A.; Gervasio, F. L.; Micheletti, C.; Parrinello, M. Efficient reconstruction of complex free energy landscapes by multiple walkers metadynamics. J. Phys. Chem. B 2006, 110, 3533–3539. (65) Altis, A.; Nguyen, P. H.; Hegger, R.; Stock, G. Dihedral angle principal component analysis of molecular dynamics simulations. J. Chem. Phys. 2007, 126, 244111. (66) Altis, A.; Otten, M.; Nguyen, P. H.; Hegger, R.; Stock, G. Construction of the free energy landscape of biomolecules via dihedral angle principal component analysis. J. Chem. Phys. 2008, 128, 245102. (67) Koukos, P. I.; Glykos, N. M. Grcarma: A fully automated task-oriented interface for the analysis of molecular dynamics trajectories. J. Comput. Chem. 2013, 34, 2310–2312.



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

TOC figure


Page 44 of 44

Molecular Dynamics Analysis of Binding of Kinase Inhibitors to WT

Recommend Documents