Title: Molecular Dynamics
1Molecular 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
2Protein 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) -
3Myoglobin
No pathway for O2 ? heme!
4Protein 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.
5Function from static structure?
6Dynamics necessary for function
7Snapshots
8Static and/or average structures may not be
representative of conformations critical to
function
9Theory
- 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
10Why 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
11Molecular Dynamics
- Potential function for MD1,2 sum of following
terms - U Bond Angle Dihedral van der Waals
Electrostatic
- Levitt M. Hirshberg M. Sharon R. Daggett V. Comp.
Phys. Comm. (1995) 91 215-231 - Levitt M. et al. J. Phys. Chem. B (1997) 101
5051-5061
12Molecular Dynamics
- Potential function for MD
- U Bond Angle Dihedral van der Waals
Electrostatic
13Molecular Dynamics
- Potential function for MD
- U Bond Angle Dihedral van der Waals
Electrostatic
b0
14Molecular Dynamics
- Potential function for MD
- U Bond Angle Dihedral van der Waals
Electrostatic
?0
15Molecular Dynamics
- Potential function for MD
- U Bond Angle Dihedral van der Waals
Electrostatic
F0
16Molecular Dynamics
- Potential function for MD
- U Bond Angle Dihedral van der Waals
Electrostatic
17Molecular 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
18Molecular Dynamics
- Non-bonded components of potential function
- Unb van der Waals Electrostatic
19Molecular Dynamics
- Non-bonded components of potential function
- Unb van der Waals Electrostatic
20Molecular Dynamics
- Non-bonded components of potential function
-
21Molecular Dynamics
- Non-bonded components of potential function
22Molecular 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!
23What 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.
24Molecular Dynamics
- Time dependent integration of classical equations
of motion
25Molecular Dynamics
- Time dependent integration
26Molecular Dynamics
- Time dependent integration
27Molecular Dynamics
- Time dependent integration
28Molecular Dynamics
- Time dependent integration
29Molecular Dynamics
- Time dependent integration
30Molecular Dynamics
- Time dependent integration
31Molecular 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
32Molecular Dynamics
Actual integration the equations of motion
Conserves energy Smooth, robust
33Molecular Dynamics
Determination of temperature
34Methods
- 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
35Molecular 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
36Methods
- 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
37Implementation
- 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
38Implementation
- 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
40Molecular Dynamics
- MD provides atomic resolution of native dynamics
PDB ID 3chy, E. coli CheY 1.66 Å X-ray
crystallography
41Molecular Dynamics
- MD provides atomic resolution of native dynamics
PDB ID 3chy, E. coli CheY 1.66 Å X-ray
crystallography
42Molecular Dynamics
- MD provides atomic resolution of native dynamics
3chy, hydrogens added
43Molecular Dynamics
- MD provides atomic resolution of native dynamics
3chy, waters added (i.e. solvated)
44Molecular Dynamics
- MD provides atomic resolution of native dynamics
3chy, waters and hydrogens hidden
45Molecular Dynamics
- MD provides atomic resolution of native dynamics
native state simulation of 3chy at 298 Kelvin,
waters and hydrogens hidden
46Molecular Dynamics
- MD provides atomic resolution of native dynamics
native state simulation of 3chy at 298 Kelvin,
waters and hydrogens hidden
47Molecular Dynamics
- MD provides atomic resolution of folding /
unfolding
unfolding simulation (reversed) of 3chy at 498
Kelvin, waters hydrogens hidden
48(No Transcript)
49Average may not be representive
50Dynamic cleft discovered through MD
Cytochrome b5
Storch et al., Biochem, 1995, 1999a,b, 2000
51Construction of mutants to test whether
cleft forms
R47
S18
Storch et al., Biochem, 1999 Bill Atkins,
Patricia Campbell
52Construction of cyt c cyt b5 complexes
53Changes in cyt b5 upon binding cyt c
Predicted binding surface
Change in chemical shift
Hom et al., Biochem, 2000
54Cleft allows for electron transfer through the
protein in channel lined with aromatics
Nonpolar
Polar
55Validation
- Validation, how do you know if a simulation is
correct? - How do you know it is done?
56Starting 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
57Native 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)
58Structural Changes to Native State During MD
Turn
Active Site Loop
N-terminus
Xtal Structure 50 ns MD
Crystal Structure
After 50 ns of MD
59Chymotrypsin Inhibitor 2
- Expt Shaw et al (1995) Biochem 342225
- MD Li Daggett (1995) Prot. Eng. 8(11)
no mobility
high
N15H
60Temperature (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
61Pushing to high temperature
- Taking excursions farther from the native state,
will the force fields and methods hold? - Thermal unfolding of proteins
62Thermal Denaturation of CI2 at 498 K
D
Ca RMSD (Å)
TS?
N
Time (ns)
Li Daggett, 1994, PNAS JMB, 1996
63Identifying 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
64Projection 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
65Main-chain Fold Preserved in Transition State
b2
a
b1
b3
Crystal Structure
Average TS1 Structure 4 Å, 43 native H-Bonds
66Packing is Disrupted in Transition State
WT TS1 32 SASA
Crystal Structure
67Structure 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,
68Calculation 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
69Comparison 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
70Overall 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
71Conformational Heterogeneity of TS
Rotate 90
to right
Crystal Structure TS1-4, 498 K TS5-9, DT
ltRMSDgtXTAL?ltRMSDgtMD ? 4.5 Å
72Free 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
73FEP 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
74Designing 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
75Removal 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.
76Thermal Denaturation of CI2 at 498 K
D 40,000 structures
Ca RMSD (Å)
Time (ns)
77The 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
78Summary 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.
79En-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)
81Development of information-rich property space
Low information content property Main-chain
non-polar SASA No discrimination between native
non-native states
Native
Nonnative
82Development 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.
83Development 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