Molecular Dynamics

1
Molecular Dynamics
...everything that living things do can be
understood in terms of the jigglings and
wigglings of atoms. Richard Feynman
Valerie Daggett Bioengineering Department Universi
ty of Washington
2
Protein Dynamics

• Proteins are not static
• Motion is an incontrovertible consequence of
  existing _at_ room temperature (or any T gt 0 K)
• Kinetic energy per atom is 1 kcal/mole _at_ 298K
  (25C) ? several Å/ps
• Motion recognized to be important early on.
  Kendrew (1950s) solved crystal structure of
  myoglobin (Perutz, phasing)

3
Myoglobin
No pathway for O2 ? heme!
4
Protein Dynamics

• Kendrew Perhaps the most remarkable features
• of the molecule are its complexity and its lack
  of
• symmetry. The arrangement seems to be almost
• totally lacking in the kind of regularities
  which one
• instinctively anticipates, and it is more
  complicated
• than has been predicted by any theory of protein
  structure.
• Situation gets worse when you consider dynamics.
• But proteins are dynamic and dynamic behavior
  critical for function. So, static, average
  structures are only part of the story.

5
Function from static structure?
6
Dynamics necessary for function
7
Snapshots
8
Static and/or average structures may not be
representative of conformations critical to
function
9
Theory

• Experiment clearly demonstrates that proteins are
  mobile, but no single experiment or combination
  of experiments can provide an all-inclusive view
  of the dynamic behavior of all atoms in a
  protein.
• Computer simulations can however
• ea. atom as a function of time

10
Why Molecular Dynamics?

• Most realistic simulation method available
• Can provide structural and dynamic information
  unobtainable by experiment, but is experimentally
  testable
• Native and nonnative interactions apparent
• But
• Sampling is limited, the goal is to sample
  experimentally relevant regions of conformational
  space, not all of conformational space

11
Molecular Dynamics

• Potential function for MD1,2 sum of following
  terms
• U Bond Angle Dihedral van der Waals
  Electrostatic

1. Levitt M. Hirshberg M. Sharon R. Daggett V. Comp.
   Phys. Comm. (1995) 91 215-231
2. Levitt M. et al. J. Phys. Chem. B (1997) 101
   5051-5061

12
Molecular Dynamics

• Potential function for MD
• U Bond Angle Dihedral van der Waals
  Electrostatic

13
Molecular Dynamics

• Potential function for MD
• U Bond Angle Dihedral van der Waals
  Electrostatic

b0
14
Molecular Dynamics

• Potential function for MD
• U Bond Angle Dihedral van der Waals
  Electrostatic

?0
15
Molecular Dynamics

• Potential function for MD
• U Bond Angle Dihedral van der Waals
  Electrostatic

F0
16
Molecular Dynamics

• Potential function for MD
• U Bond Angle Dihedral van der Waals
  Electrostatic

17
Molecular Dynamics

• Non-bonded components of potential function
• Unb van der Waals Electrostatic
• To a large degree, protein structure is dependent
  on non-bonded atomic interactions

18
Molecular Dynamics

• Non-bonded components of potential function
• Unb van der Waals Electrostatic

19
Molecular Dynamics

• Non-bonded components of potential function
• Unb van der Waals Electrostatic

20
Molecular Dynamics

• Non-bonded components of potential function

-
21
Molecular Dynamics

• Non-bonded components of potential function

22
Molecular Dynamics

• Non-bonded components of potential function

NOTE Sum over all pairs of N atoms, or
pairs
N is often between 5x105 to 5x106 For 5x105 that
is 1.25x1011 pairs THAT IS A LOT OF POSSIBLE
PAIRS!
23
What can you do with a force field? Generation
of experimental structures Refinement of
experimental structures Monte Carlo Scoring
functions Energy minimization Analysis Perform MD
simulations etc.
24
Molecular Dynamics

• Time dependent integration of classical equations
  of motion

25
Molecular Dynamics

• Time dependent integration

26
Molecular Dynamics

• Time dependent integration

27
Molecular Dynamics

• Time dependent integration

28
Molecular Dynamics

• Time dependent integration

29
Molecular Dynamics

• Time dependent integration

30
Molecular Dynamics

• Time dependent integration

31
Molecular Dynamics

• Time dependent integration

Evaluate forces and perform integration for every
atom Each picosecond of simulation time requires
500 iterations of cycle E.g. w/ 50,000 atoms,
each ps (10-12 s) involves 25,000,000 evaluations
32
Molecular Dynamics
Actual integration the equations of motion
Conserves energy Smooth, robust
33
Molecular Dynamics
Determination of temperature
34
Methods

• Molecular dynamics (MD)
• Brooks-Beeman integration algorithm
• Microcanonical ensemble (NVE)
• Number of atoms, box volume, energy are
  conserved
• Energy conservation is naturally satisfied with
  classical equations of motion
• Energy conservation is an inherent check on the
  implementation
• Free from coupling the microscopic system to
  macroscopic variables as do NVT and NPT

35
Molecular dynamics

• Microcanonical ensemble, all atoms, solvent,
• fully flexible molecules, continuous
  trajectories,
• no restraints/biases
• predictive MD---expt to check
• No Ewald --- artificial periodicity, altered
  conformational and dynamical properties
• No fictitious bonds between H atoms of water
• No Shake
• Correct masses
• Good simple, flexible water correct D and RDF

36
Methods

• Molecular dynamics (MD)
• Temperature in NVE
• Mean T over hundreds of steps
• Energy drift in ilmm is primarily kinetic
  resulting from numerical round-off
• Over thousands of steps mean T can be monitored
  for energy conservation
• Velocity rescales once per 10ns

37
Implementation

• Written in C
• Ubiquitous, standardized, optimized language
• 64 bit math
• Software design
• Kernel
• Compiles users molecular mechanics programs
• Schedules execution across processor and machines
• Modules, e.g.
• Energy minimization
• Molecular Dynamics
• Monte Carlo
• Analysis
• REMD
• RDCs
• others

38
Implementation

• Dual mode parallelization
• Standardized tools available on modern platforms
• POSIX threads
• Distribute computations across multiple CPUs in a
  single computer
• Message Passing Interface (MPI)
• Distribute computations across multiple computers
  on a high speed network
• Benefit is scalability

39

• State of the Art MD
• What can be done with PCs?
• Environment
• Possible to characterize solvent-dependent
• conformational behavior
• Proteins in membranes
• Size
• 500 residues (more possible if willing to
  dedicate resources
• to it, our record is 2519 residues in
  solvated membrane)
• Timescale
• Multiple 20-100 ns simulations fairly routine
  for proteins
• ms possible if willing to dedicate resources to
  it

40
Molecular Dynamics

• MD provides atomic resolution of native dynamics

PDB ID 3chy, E. coli CheY 1.66 Å X-ray
crystallography
41
Molecular Dynamics

• MD provides atomic resolution of native dynamics

PDB ID 3chy, E. coli CheY 1.66 Å X-ray
crystallography
42
Molecular Dynamics

• MD provides atomic resolution of native dynamics

3chy, hydrogens added
43
Molecular Dynamics

• MD provides atomic resolution of native dynamics

3chy, waters added (i.e. solvated)
44
Molecular Dynamics

• MD provides atomic resolution of native dynamics

3chy, waters and hydrogens hidden
45
Molecular Dynamics

• MD provides atomic resolution of native dynamics

native state simulation of 3chy at 298 Kelvin,
waters and hydrogens hidden
46
Molecular Dynamics

• MD provides atomic resolution of native dynamics

native state simulation of 3chy at 298 Kelvin,
waters and hydrogens hidden
47
Molecular Dynamics

• MD provides atomic resolution of folding /
  unfolding

unfolding simulation (reversed) of 3chy at 498
Kelvin, waters hydrogens hidden
48
(No Transcript)
49
Average may not be representive
50
Dynamic cleft discovered through MD
Cytochrome b5
Storch et al., Biochem, 1995, 1999a,b, 2000
51
Construction of mutants to test whether
cleft forms
R47
S18
Storch et al., Biochem, 1999 Bill Atkins,
Patricia Campbell
52
Construction of cyt c cyt b5 complexes
53
Changes in cyt b5 upon binding cyt c
Predicted binding surface
Change in chemical shift
Hom et al., Biochem, 2000
54
Cleft allows for electron transfer through the
protein in channel lined with aromatics
Nonpolar
Polar
55
Validation

• Validation, how do you know if a simulation is
  correct?
• How do you know it is done?

56
Starting a Molecular Dynamics Simulation
Heat to desired temperature and allow motion to
evolve over time
Solvate with water or other solvent 8 - 14
Å from protein
T 298 K r 0.997 gm/ml T 498 K r 0.829
gm/ml
Crystal or NMR Structure
All atoms present Fully flexible water NVE
57
Native Dynamics at 25 ºC
3
All Ca atoms
ltRMSDgt 1.7 Å
2.5
2
1.5
NOEs reproduced
Ca RMSD (Å)
1
Active site loop and N-terminus removed
0.5
ltRMSDgt 0.7 Å
0
Crystal structure
5
15
40
35
30
25
20
50
45
10
0
Time (ns)
58
Structural Changes to Native State During MD
Turn
Active Site Loop
N-terminus
Xtal Structure 50 ns MD
Crystal Structure
After 50 ns of MD
59
Chymotrypsin Inhibitor 2

• Expt Shaw et al (1995) Biochem 342225
• MD Li Daggett (1995) Prot. Eng. 8(11)

no mobility
high
N15H
60
Temperature (K) Cutoff range (Å) Timea (ns) NOEs satisfiedb ( of 603)
Xtal 91.00
298 8 100 98.18
298 8 50 98.51
298 8 50 98.00
298 8 20 97.18
298 8 20 97.84
298 8 20 98.18
298 10 100 98.84
298 10 50 97.68
298 10 50 97.51
298 10 20 97.84
298 10 20 97.54
298 10 20 97.35
310 8 100 98.51
310 8 50 98.01
310 8 50 98.51
310 10 100 98.68
310 10 50 98.01
310 10 50 98.51
323 8 100 98.34
323 8 50 97.35
323 10 100 98.34
323 10 50 98.01
No restraints
22 simulations gt1.2 ms
NOEs courtesy of Stefan Freund Trevor
Rutherford, ARF
61
Pushing to high temperature

• Taking excursions farther from the native state,
  will the force fields and methods hold?
• Thermal unfolding of proteins

62
Thermal Denaturation of CI2 at 498 K
D
Ca RMSD (Å)
TS?
N
Time (ns)
Li Daggett, 1994, PNAS JMB, 1996
63
Identifying Transition States in MD Trajectories

• TS not localized to single bond, distributed and
  ensemble
• From MD cannot calculate DG along reaction
  coordinate
• Structure-based definition of TS
• Kinetically, protein will not succeed in every
  attempt to cross TS but will change rapidly
  afterwards
• A process with a large change in energy but small
  change in entropy ? large change in free energy

H
TS
G
N
-TS
D
Reaction Coordinate
64
Projection of Trajectory in RMSD Space
D
TS
N
Calculate the RMSD between all structures---15,000
x 15,000 dimensional space Reduce to
3-dimensions, distance between points ? RMSD
between structures Clusters indicate similar
conformations, conformational states
Li Daggett, 1994, 1996
65
Main-chain Fold Preserved in Transition State
b2
a
b1
b3
Crystal Structure
Average TS1 Structure 4 Å, 43 native H-Bonds
66
Packing is Disrupted in Transition State
WT TS1 32 SASA
Crystal Structure
67
Structure of TS from Experiment
N
TS
D
TS
TS
DDGTS-D
DDGTS-D
D
D
N
DGN-D
N
DGN-D
DDGN-D
DDGN-D
F DDGTS-D / DDGN-D 0
F DDGTS-D / DDGN-D 1
F 1 ? site of mutation native-like in TS F 0
? site of mutation unfolded in TS Fractional F
values ? partial structure in TS
Matouschek, Kellis, Serrano, Fersht, 1989,
Nature, Fersht, Leatherbarrow, Wells, 1986,
Nature 1987 Biochem,
68
Calculation of S Values for Comparison with
Experimental FF Values
For each residue calculate
S structure index (S2º ) (S3º )
native secondary structure (f,j)
tertiary structure, contacts
Daggett Li, 1994, PNAS Daggett, Li, Itzhaki,
Otzen Fersht, 1996, JMB
69
Comparison of Calculated S Values and
Experimental F Values
a
b1
b2
b3
Phi or S Value
S Value
Residue Number
Phi Value
Otzen et al., 1994, PNAS Itzhaki et al., 1995,
JMB Li Daggett, 1994,1996, JMB Daggett et
al., 1996, JMB
70
Overall TS Structure and Unfolding Pathway are
Independent of Temperature
TS1
373 K
(0.225 ns)
(21 ns)
398 K
TS2
(8.26 ns)
(0.335 ns)
448 K
TS3
(1.44 ns)
(0.1 ns)
473 K
TS4
(0.57 ns)
(0.07 ns)
498 K
(0.3 ns)
498 K
DT
71
Conformational Heterogeneity of TS
Rotate 90
to right
Crystal Structure TS1-4, 498 K TS5-9, DT
ltRMSDgtXTAL?ltRMSDgtMD ? 4.5 Å
72
Free Energy Calculations for Direct Determination
of F
TS
N
D
DGN?D
DGN?TS
DGTS?D
WT
N
TS
D
DGN
DGTS
DGD
DG'TS?D
DG'N?TS
Mutant
N'
TS'
D'
DG'N?D
73
FEP Calculations for Hydrophobic Core Mutants
R 0.85 R 0.91 (no V47A)
(kcal/mole)
Extended peptide NOT a good model of D
Pan Daggett, Biochem, 2001
74
Designing Faster Folding Forms of CI2 Based on
MD-Generated TS Models
F50
F50
F48
R62
E7
R62
R48
E7
K2
K2
D23
A23
DA 23 TS
WT TS
RF48 TS
Ladurner, Itzhaki, Daggett, Fersht, 1998, PNAS
75
Removal of Unfavorable Interactions Identified
in TS Models Accelerates Folding
Removal of charge repulsion and improvement of
packing in the TS yields fastest-folding form of
CI2.
76
Thermal Denaturation of CI2 at 498 K
D 40,000 structures
Ca RMSD (Å)
Time (ns)
77
The Denatured State of CI2
hydrophobic clustering
Distances in N W5-V14 15 Å I30-Y42 13 Å P33-I37
11 Å
P33
I37
I30
Experimental Results
lt3JNH-CaHgtexpt 7.2 Hz
Y42
lt3JNH-CaHgtMD 7.0 Hz
a
L49
(res. 17-21)
I57
V14
Dd V19-L21, I30-T36 Nearly random coil
W5
Kazmirski et al., PNAS, 2001
78
Summary of CI2 Simulations

• N is well behaved and in good agreement with
  experiment.
• TS is an expanded version of N with disrupted
  core and loops and frayed secondary structure.
• Validity of MD-generated TS models tested through
  indirect comparison with experimental F values,
  direct comparison of DGs, behavior when T is
  quenched, and design of faster folding mutants.
• WT D is very disrupted with only minor amounts of
  hydrophobic clustering and fluctuating helical
  structure. Nearly random coil.

79
En-HD unfolds at 348 373 K on the same
timescale by simulation and experiment
Time to reach TS in MD simulation
kunf
47,000 s-1
10 C
But, it is not enough to get the
timescale right, must get pathway too!
kf
Mayor et al., Nature 2003
80
(No Transcript)
81
Development of information-rich property space
Low information content property Main-chain
non-polar SASA No discrimination between native
non-native states
Native
Nonnative
82
Development of information-rich property space
High information content property CONGENEAL
structural dissimilarity score1 Excellent
discrimination between native non-native states
Native
Nonnative
1. Yee and Dill, Protein Science, 1993.
83
Development of a reaction coordinate from
property space

• Foldedness ? location along folding reaction
  coordinate

• Mean distance in PS (32 -gt 10 properties) for a
  given conformation to the folded or native state
  cluster is acceptable reaction coordinate
• Value increases with distance from native cluster
• Native cluster is bounded