Multiscale Modeling of Phosphorylation and Inhibition of the Epidermal Growth Factor Receptor Tyrosine Kinase: Linking Somatic Mutations to Differential Signaling

Multiscale Modeling of Phosphorylation and Inhibition of the Epidermal Growth Factor Receptor Tyrosine Kinase: Linking Somatic Mutations to Differential Signaling

1. University of Pennsylvania Department of Bioengineering
Multiscale Modeling of Phosphorylation and
Inhibition of the Epidermal Growth Factor
Receptor Tyrosine Kinase: Linking Somatic
Mutations to Differential Signaling
Yingting Liu
Advisor: Dr. Ravi Radhakrishnan
Department of Bioengineering
University of Pennsylvania

2. University of Pennsylvania Department of Bioengineering
Outline
Backgrounds
Hypothesis and Specific aims
Experimental design and preliminary results

3. University of Pennsylvania Department of Bioengineering
ErbB Family Receptors and the Signaling Pathways
Yarden and Sliwkowski, nature reviews, 2001

4. University of Pennsylvania Department of Bioengineering
Tyrosine Phosphorylation and Receptor Inhibition
Zhang and Kuriyan,Cell, 2006

5. University of Pennsylvania Department of Bioengineering
EGFR Kinase Domain Mutations
Choi and Lemmon, Oncogene, 2007 Zhang and Kuriyan,Cell, 2006Carey and Sliwkowski, Cancer Res, 2006

6. University of Pennsylvania Department of Bioengineering
Hypothesis and Methods
WehypothesizethatmutantsintheEGFR
kinasedomainwillalterthekinase-inhibitor,
kinase-substrateinteractions,andthe
catalyticreactionefficiencyoftheturn-over
ofdifferentEGFRsubstratesbyaffecting
thepropertiesofEGFRTKactivesite,
thereforeleadtodifferentialcharacteristics
inthedownstreamsignalinginpathways
mediatedbyEGFR.
Weproposetoemploymultiscale
computationalmethodsbasedonmolecular
docking,moleculardynamics(MD),and
quantummechanicsmolecularmechanics
(QM/MM) simulationstotestthis
hypothesis.MD simulation for protein kinase
Multiple conformation molecular docking
MD simulation for complex
Structural and energetic analysis
Inhibition
QM/MM calculation on catalysis
MD simulation for
EGFRTK-ATP-MG-Peptide
complex
MD simulation for protein kinase
Multiple conformation molecular docking
MD simulation for complex
Structural and energetic analysis
Inhibition
QM/MM calculation on catalysis
MD simulation for
EGFRTK-ATP-MG-Peptide
complex

7. University of Pennsylvania Department of Bioengineering
Specific Aims
Aim1.Developingempiricalforce-fieldparametersforsmallmolecule
inhibitorsforuseinin-silicodockingandmoleculardynamicssimulations.
Aim2.Exploringtheconformationalandfreeenergylandscapefor
wildtypeandL834RmutantEGFRkinasecomplexedwithsmallmolecule
inhibitorsandpeptidesubstrates.
Aim3.ModelingthecatalyticmechanismandactivityoftheEGFR
tyrosinekinase.

8. University of Pennsylvania Department of Bioengineering
Specific Aims
Aim1.Developingempiricalforce-fieldparametersforsmallmolecule
inhibitorsforuseinin-silicodockingandmoleculardynamics
simulations.
Aim2.Exploringtheconformationalandfreeenergylandscapefor
wildtypeandL834RmutantEGFRkinasecomplexedwithsmall
moleculeinhibitorsandpeptidesubstrates.
Aim3.ModelingthecatalyticmechanismandactivityoftheEGFR
tyrosinekinase.

9. University of Pennsylvania Department of Bioengineering
MD Simulation and CHARMM Potential EnergyV 2 2 2
( ) ( - ) ( - ) ( - ) (1 cos( - ))
0 0 0
12 6
min min2
( - ) -
0
V K b b K S S K K n
b UB
bonds UB angles dihedrals
R R q q
ij ij i j
K j j
imp ij r r er
impropers nonbond
ij ij ij
           

   

   

   

   

   

r
Molecular Dynamic (MD) simulations:
CHARMM potential energy:
Essential part is the potential energy function.

10. University of Pennsylvania Department of Bioengineering
Erlotinib Parameterization (1)Q
2AN
C
AQ
1
N
AQ
1
CAQ
3
CAQ
3
Q
2AC
C
A
CAQ
4
CAQ
4
CA
S
O
O
S
C
T
2
T
2C
O
S
T
2C
T
2
C
S
O
T
3C
T
3C
CA
CA
CA
CA
CA
CA
CC
3
CC
3
HA
H
PH
PH
HP
HP
HP
PH
HA
HA
HA
HA
HA
HA
HA
HA
HA
HA
AH
HP
HA
AQ
1
N
Define new atom types and initiate the parameter set.
Optimize the structure using ab-initio methods and obtain
equilibrium constants.
Obtain partial charges of each atom using CHELPG
(CHarges from ELectrostatic Potentials using a Grid based method).
Get Van der Waals constants ( and ) from existing
CHARMM parameters.
Guess the force field constants based on those assigned for
similar structure in existing CHARMM parameters.minij
R ij
 2 2 2
( ) ( - ) ( - ) ( - ) (1 cos( - ))
0 0 0
12 6
min min2
( - ) -
0
V K b b K S S K K n
b UB
bonds UB angles dihedrals
R R q q
ij ij i j
K j j
imp ij r r er
impropers nonbond
ij ij ij
           

   

   

   

   

   

r

11. University of Pennsylvania Department of Bioengineering
Erlotinib Parameterization (2)N
3
C
19
N
2
C
18
C
7
C
6
C
17
C
13
C
9
C
8
N
1
N
1H
H
8
H
19
H
17
H
22
O
2
H
21
H
11
O
1
H
12
O
3
H
33
H
32
Refine Partial charges manually.Q
2AN
C
AQ
1
N
AQ
1
CAQ
3
CAQ
3
Q
2AC
C
A
CAQ
4
CAQ
4
CA
S
O
O
S
C
T
2
T
2C
O
S
T
2C
T
2
C
S
O
T
3C
T
3C
CA
CA
CA
CA
CA
CA
CC
3
CC
3
HA
H
PH
PH
HP
HP
HP
PH
HA
HA
HA
HA
HA
HA
HA
HA
HA
HA
AH
HP
HA
AQ
1
N

12. University of Pennsylvania Department of Bioengineering
Erlotinib Parameterization (3)
Refine dihedral parameters to reproduce ab initio dihedral energy surface.
Using genetic algorithm to automatically minimize the merit function: Q
2AN
C
AQ
1
N
AQ
1
CAQ
3
CAQ
3
Q
2AC
C
A
CAQ
4
CAQ
4
CA
S
O
O
S
C
T
2
T
2C
O
S
T
2C
T
2
C
S
O
T
3C
T
3C
CA
CA
CA
CA
CA
CA
CC
3
CC
3
HA
H
PH
PH
HP
HP
HP
PH
HA
HA
HA
HA
HA
HA
HA
HA
HA
HA
AH
HP
HA
AQ
1
N 2
1
()
NGRID
CG
ii
i
DD

 0 100 200 300 400
0
0.5
1
1.5
Dihedral (degree)
Energy (Kcal/mol)
Dihedral potential energy surface C
i
D G
i
D
NGRID is the number of potential values
calculated in the surface. and
are potential values from CHARMM and
GAUSSIAN.

13. University of Pennsylvania Department of Bioengineering
Erlotinib Parameterization (4)Q
2AN
C
AQ
1
N
AQ
1
CAQ
3
CAQ
3
Q
2AC
C
A
CAQ
4
CAQ
4
CA
S
O
O
S
C
T
2
T
2C
O
S
T
2C
T
2
C
S
O
T
3C
T
3C
CA
CA
CA
CA
CA
CA
CC
3
CC
3
HA
H
PH
PH
HP
HP
HP
PH
HA
HA
HA
HA
HA
HA
HA
HA
HA
HA
AH
HP
HA
AQ
1
N
Refine force constants to reproduce vibrational eigenvalues and eigenvectors.
Using genetic algorithm to automatically minimum the merit function: 36
2
max
1
()
36
N
CG
i i j
i
N
  






 1
max ( )
i cG
j i j




Vaiana, Computer Physics Communications, 2005. , : the ith frequency and eigenvector from CHARMM and GAUSSIAN
Project into { }, and find the index jmax which maxmum .
In the ideal case, and
max
;,

  
GG
ij
CC
ii
c G c G
i j i j
C G c G
i j i j ij

   
    


14. University of Pennsylvania Department of Bioengineering
Erlotinib Parameterization (5)
------Preliminary
results
Water interactionsInteraction Energies
(Kcal/mol)
Distance (Å)
GAUSSIAN CHARMMGAUSSIAN CHARMM
N2…HOH -6.69 -6.61 2.13 1.91
N3…HOH_2 -5.33 -5.3 2.32 2.01
N1H…OHH_2 -6.52 -6.52 2.4415 2.63
Dipole moment
(Debye)
GAUSSIAN CHARMM
4.868 5.07
Table 1 Water-mediated interactions and dipole moment for erlotinib. The
ab-initio interaction energies are scaled by 1.16, and the distances should
offset by –0.1 to –0.2 A. Experimental dipole moments are typically ~10 to
20% larger than HF/6-31G*.

15. University of Pennsylvania Department of Bioengineering
Erlotinib Parameterization (6)
------Preliminary
results0 50 100 150 200 250 300 350
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
Dihedral (degree)
Energy (kcal/mol)
Frequency matching Potential surface fitting
Genetic algorithm efficiency

16. University of Pennsylvania Department of Bioengineering
Specific Aims
Aim1.Developingempiricalforce-fieldparametersforsmallmolecule
inhibitorsforuseinin-silicodockingandmoleculardynamicssimulations.
Aim2.Exploringtheconformationalandfreeenergylandscapefor
wildtypeandL834RmutantEGFRkinasecomplexedwithsmallmolecule
inhibitorsandpeptidesubstrates.
Aim3.ModelingthecatalyticmechanismandactivityoftheEGFR
tyrosinekinase.

17. University of Pennsylvania Department of Bioengineering
Methods: MD simulations
Solvated model for MD simulation of EGFRTK.
(Iceblue: sodium; yellow: chlorine; orange: protein; tan: water).
MolecularDynamic(MD)protocol:
•Prepareproteinconformationbased
onavailablecrystalstructureor
homologies.
•Solvatetheproteinandneutralize
thesystemsbyplacingionsrandomly.
•Minimizethesolvatedmodels
•Heatthesystemto300K
•Equilibrateatconstanttemperature
andconstantpressure(300Kand1
atm)for200pstostablethesystem.
•Runproductivetrajectory.

18. University of Pennsylvania Department of Bioengineering
Methods: Multiple-Conformation Molecular Docking
•Theideaofmoleculardocking:togenerateacomprehensivesetof
conformationsofthereceptor-ligandcomplexandthentorankthem
accordingtotheirstability.
•Singleconformationdocking:Ligandisflexibility,whilereceptors
areusuallytreatedasrigidduringdocking.
•Multiple-conformationdocking:Anensembleof100snapshotsof
theproteiniscollectedfromtheequilibratedtrajectorytoperform
moleculardocking.Thegeneratedligandconformationsareclustered
basedontherelativeRMSDandanalyzedtoexplorethe
conformationalandfreeenergylandscapeoftheinteractionbetween
proteinkinaseandtheligands.
Themultiple-conformationdockingjobsaresubmittedinparallelso
thattheywillrunsimultaneouslyandthenclusterthegenerated
conformationsuponcompletionofthedockingrunsusingFortran90
program.

19. University of Pennsylvania Department of Bioengineering
Methods: Binding Free Energy Calculation  
22
12 6 12 6
( / 2 )
()
()

ij
ij ij ij ij i j
vdW hbond elec
ij ij ij ij ijij ij ij ij
r
tor tor sol i j j i
ij
A B C D q q
G G G E t G
rrR R R R
G N G SV S V e



   
           
   
   
     
  

•MolecularMechanicsPoisson-BoltzmannSurfaceArea(MMPBSA): - - ;
- ;
;
.
complex receptor ligand
MM PBSA MM
MM bond angle tors elec vdw
PBSA solvation PB SA
G G G G
G E G TS
E E E E E E
G G G G


    
  
Electrostatic solvation energy: Poisson-
Boltzmann equation.
Nonpolar term: depend on surface area.
Sitkoff and Honig 1993
•AUTODOCK:

20. University of Pennsylvania Department of Bioengineering
Kinase-Inhibitor Interactions
------Proposed model
•Motivation:Similarbindingconformationspresentedincrystalstructuresbut
remarkablyincreasethebindingaffinitiesinL834Rmutantsystems.
---erlotinib(CareyandSliwkowski,2006),gefitinibandAEE788(YunandEck
2007)
•SpecificofAim:usingmultiple-conformationmoleculardockingtoobtainsix
toprankedcomplexconformationsbasedontheapproximatefreeenergyfrom
AUTODOCK andthenperformMDbasedstructuralandenergeticanalysis
(MMPBSA)foreachconformations.Amongthesix,threeconformationswillbe
highlightedforanalysisbasedonthemoreaccuratebindingfreeenergy.
•Possiblereasonstotest:uniqueinteractionsbetweenL834Rmutantkinaseand
inhibitors,subtleconformationaldifferences,whichishardtobecapturedby
crystallographicmethods,effectofsolvation,…

21. University of Pennsylvania Department of Bioengineering
Kinase-Inhibitor Interactions
------Preliminary results and future work
WT L858R
Crystal
conf.
Lowest energy conf.
Top ranked Erlotinib conformations in EGFR wildtype and mutant system.
Use MD simulations to refine these structures with explicit solvent and resort the
structures using MMPBSA methods.
Perform structural analysis for the refined conformations to explore the effect of
mutations on kinase-inhibitor interaction.

22. University of Pennsylvania Department of Bioengineering
Kinase-Substrate interactions
------Proposed model
•SpecificofAim:performthemultiple-conformationmoleculardocking
protocolfollowedbytheMDbasedstructuralanalysisandfreeenergy
calculationtopredictthebestbindingmodesandobtainthecorresponding
bindingaffinities,whichcanbecorrelatedtoKmvaluesforeachsubstrate.
•Motivation:topredictthebinding
modesfordifferentsubstratesandtest
theeffectofmutationonkinase-
substrateinteraction.
•Substrates:Fourseven-residue
sequencesderivedfromtheC-terminal
tailofEGFRTK(Y1068,Y1173,Y992
andY1045).

23. University of Pennsylvania Department of Bioengineering
Kinase-Substrate interactions
------Preliminary results and future work
L858R unphosphorylated EGFR
Binding with 10687
2GS6

24. University of Pennsylvania Department of Bioengineering
Kinase-Subtrate interactions
------Preliminary results and future work
Substrate
s
Approximate Binding
Energy(Kcal/mol)
Y1068 Y1173 Y992
Wildtype-5.42 -4.69 -4.7
L834R
mutant
-5.93 -3.78 -5.91
Liu, Purvis and Radhakrishnan,2007

25. University of Pennsylvania Department of Bioengineering
Kinase-Substrate interactions
------Preliminary results and future work
FreeenergycontributionsofEGFRTK-peptide
(VPEYINQ)bindingfromMMPBSA calculation.
(Kcal/mol)
Internal energy -139.7
Polar solvation 140.5
onpolar solvation -6.4
Total binding free energy -5.6
Future work: Use MD simulations to refine these structures with explicit
solvent and recalculate the binding free energy using MMPBSA methods.

26. University of Pennsylvania Department of Bioengineering
Specific Aims
Aim1.Developingempiricalforce-fieldparametersforsmallmolecule
inhibitorsforuseinin-silicodockingandmoleculardynamicssimulations.
Aim2.Exploringtheconformationalandfreeenergylandscapefor
wildtypeandL834RmutantEGFRkinasecomplexedwithsmallmolecule
inhibitorsandpeptidesubstrates.
Aim3.ModelingthecatalyticmechanismandactivityoftheEGFR
tyrosinekinase.

27. University of Pennsylvania Department of Bioengineering
Catalytic Mechanism
•In principle, the reaction mechanism can be either an associative or dissociative
pathway.
•pKa and nucleophile coefficient measurements support a dissociative transition
state. (Kim and Cole, 1998)
•QM/MM studies of cAMP agree with the dissociate mechanism. (Cheng and
McCammon, 2005)

28. University of Pennsylvania Department of BioengineeringCH
2
O
P
O
-
O
O
P
O
O
O
-
P
O
-
O
-
Mg
2+
Mg
2+
O
H H2
C
O
O
O
-
CH
2
ATP
Peptide
Asp813 CH
2
O
P
O
-
O
O
P
O
O
O
-
P
O
-
O
Mg
2+
Mg
2+
O
H H2
C
O
O
-
CH
2
O
-
Asp
813
Peptide
ATP CH
2
O
P
O
-
O
O
P
O
O
O
-
P
O
-
O
Mg
2+
Mg
2+
O
H H2
C
O
O
-
CH
2
O
Asp
813
Peptide
ATP
Proposed Catalytic Mechanism for EGFRTK based on cAMP

29. University of Pennsylvania Department of Bioengineering
Prepare the Enzyme-Substrate System
Blue: 2GS6 bisubstrate;
Pink: ATP conformation in
2ITX;
Yellow: proposed peptide
conformation in aim 2;

30. University of Pennsylvania Department of Bioengineering
QM/MM Calculation
Molecular Mechanics (MM): cannot account for the covalent transformations of
chemical bonds.
Quantum Mechanics (QM): limited system size due to computational complexity.
QM/MM: Treat atoms involved in chemical reaction with QM and others MM.
QM region
MM region
Link atoms
ATP
PEPTIDE
MG

31. University of Pennsylvania Department of Bioengineering
Umbrella Sampling
•Umbrellasamplingenablesthecalculationofthepotentialofmeanforce(free
energydensity)alonganapriorichosensetofreactioncoordinates(ororder
parameters),fromwhichfreeenergychangesarecalculatedbynumerical
integration.( ) ( ) ( )u r u r W r 2
0
( ) ( )
w
W r k r r ()ur ()ur

32. University of Pennsylvania Department of Bioengineering
Free Energy Landscape Along Reaction CoordinatesCH
2
O
P
O
-
O
O
P
O
O
O
-
P
O
-
O
-
Mg
2+
Mg
2+
O
H H2
C
O
O
O
-
CH
2
ATP
Peptide
Asp813 1
r 2
r
•Umbrella sampling along two coordinates.
•25 windows are sampled as a uniform 5×5 grid along
r
1
and r
2
.
•with each window harvesting a QM/MM MD
trajectory of 2 ps.
•free energy profile as a function of the coordinate will
be calculated using the weighted histogram analysis
method (WHAM).
•Explore the effect of mutation on the reaction profile.
Gregersen and York 2003

33. University of Pennsylvania Department of Bioengineering
Summary and Significances
Effect of mutation on:
•Kinase-Inhibitor interactions.
•Kinase-Substrate interactions.
•EGFR tyrosine kinase reaction profile.
Significances:
•generatearichamountofinformationconcerningstructuralanddynamic
propertiesofthesystematatomiclevel.
•helptofurtherunderstandthemechanismofproteinkinasesinhibition
andphosphorylationandthereforeguidecancertherapyofproteinkinase
systems.

34. University of Pennsylvania Department of Bioengineering
Thanks.

35. University of Pennsylvania Department of Bioengineering
Mutations increase kinase activities
Yun et al., (Eck) Cancer Cell (2007)
Zhang et al., (Kuriyan) Cell (2006)

36. University of Pennsylvania Department of Bioengineering
Structural Studies of EGFRTK Active Site
αC-helix
peptideC-loop
ATP
GLU738
LYS721
ASP813
ASP831
MET769
G-loop
N-lobe
C-lobe
A-loop