Biophysical Chemistry G4170:

1
Biophysical Chemistry G4170 Introduction
to Molecular Dynamics
Ruhong Zhou
IBM Thomas Watson Research Center Yorktown
Heights, NY 10598
2
Molecular Simulations

• NMR or X-ray structure refinements
• Protein structure prediction
• Protein folding kinetics and mechanics
• Conformational dynamics
• Global optimization
• DNA/RNA simulations
• Membrane proteins/lipid layers simulations

3
Contents

• Overview
• Force Fields
• Solvation Models
• Molecular Dynamics
• Advanced Topics

4
III. Molecular Dynamics Monte Carlo
5
Jiggling and Wiggling

• Feynman Lectures on Physics
• Certainly no subject or field is making more
  progress on so many fronts at the present moment
  than biology, and if we were to name the most
  powerful assumption of all, which leads one on
  and on in an attempt to understand life, it is
  that all things are made of atoms, and that
  everything that living things do can be
  understood in terms of the jigglings and
  wigglings of atoms.
• -Feynman, 1963

6
Historical Perspective on MD
7
Anti-Microbial Peptide Simulation
8
Molecular Dynamics

• Solve Newtons equation for a molecular system

• Or equivalently solve the classical Hamiltonian
  equation

9
Integrator Verlet Algorithm
Start with r(t), v(t), integrate it to
r(tDt), v(tDt)
r(tDt), v(tDt)
The new position at tDt
r(t), v(t)
(1)
Similarly, the old position at t-Dt
(2)
Add (1) and (2)
(3)
Thus the velocity at t is
(4)
10
Typical MD Flowchart
Program MYMD simple MD program call
init initialization t 0 do while (t .lt.
tmax) MD loop call force (x, f,
en) calculate the force call integrate (x,
f, en) integrate equation of motion t t
delt call sample sample averages enddo
stop end
11
Cutoff Smooth Cutoff

• For pairwise non-bonded interactions, cutoffs are
  often used, particularly in old days
• Lennard-Jones potential decays as 1/r6, which is
  appropriate to use a short cutoff (9-12Å)
• However, the electrostatic interaction decays
  only as 1/r, which is a long-range interaction
  and short cutoffs might have serious artifacts
• Smooth cutoffs are often preferred to make the
  first and second derivatives of the potentials
  smooth

12
Periodic Boundary ConditionsMinimum Image
rc
Central simulation box
13
Periodic Boundary ConditionsEwald Method
Ewald Sum split into two terms, real space sum
and k-space sum
14
Ensembles

• NVE micro-canonical ensemble
• NVT canonical ensemble
• NPT grand-canonical ensemble

• Temperature control
• Berendsen thermostat (velocity rescaling)
• Andersen thermostat (velocity resampling)
• Nose-Hoover chain
• Pressure control
• Berendsen volume rescaling
• Andersen piston

15
MD Analysis

• Various trajectories of energy terms,
  temperature, and other properties
• Ramachandran plots, distance contact maps
• Secondary structures
• RMS deviations
• Radius of gyration
• Autocorrelation function
• Free energy contour maps

16
Energy Fluctuations of BPTI

• BPTI 56 residues
• solvated in explicit water 8000 atoms
• PME and cutoff for electrostatics
• NPT simulation at 300K and 1atm
• Berendsen velocity and volume scaling

17
RMSD and Rg Trajectories
18
Ramachandran plots
19
Minimum Distance Matrix
20
Autocorrelation Functions
21
Free Energy Map or PMF

• Potential of Mean Force (PMF)
• Select reaction coordinates
• Histogram configurations at the reaction
  coordinates
• Calculate the probability of particular
  coordinates along the reaction path
• Finding folding mechanism
• Select reaction coordinates can be very tricky.
  It is more of an art than a science process

22
Monte Carlo Method

• Another way to generate a canonical ensemble
  
• Monte Carlo Method
• (i). Randomly generate a new conformation. For
  proteins, usually in torsional space, such as
  rotate around a dihedral angle
• (ii). Accept or reject the new conformation based
  on the Metropolis criterion

where
23
MC vs. MD

• Only energy is needed
• Straightforward to perform NVT and NPT
• Easy to constrain some degrees of freedoms (not
  include them in trials)
• For some systems, large motions can be made (LJ
  particles) between consecutive configurations
• Hard to make trials for complex systems, such as
  proteins, since proteins move collectively
• Step size decreases with system size
• Cant easily get kinetic information

• Both energy and force are needed
• Requires temperature and pressure control for NVT
  and NPT
• Needs special techniques to constrain some
  degrees of freedoms
• The consecutive configurations are very similar
• MD can move simple and complex systems the same
  way
• Same time step can be used for small or large
  systems
• Can generate kinetic data as well as
  thermodynamic data

24
Hybrid Monte Carlo

• A subtle problem for conventional MC is that most
  of the moves might be rejected due to the protein
  atom collisions
• Protein atoms move collectively
• HMC often called bad MD but good Monte Carlo
• Three step algorithm for HMC
• Sampling velocities according to Maxwell
  distribution
• Run MD for a few steps with randomly sampled
  velocities
• Accept/reject the new conformation by Metropolis
  criterion
• HMC gives the correct canonical ensemble given
  that the MD algorithm is reversible and sympletic

25
(No Transcript)
26
Protein Folding Fast Folders
Time Scale
80s 90s 00s
00s 90s 80s
ps ns ms
ms ms sec
Folding Experiments
Folding MD Simulations

• Trp-cage, designed mini-protein (20 aa) 4ms
• b-hairpin of C-terminus of protein G (16 aa)
  6ms
• Engrailed homeodomain (En-HD) (61 aa) 27ms
• WW domains (38-44 aa) gt24ms
• Fe(II) cytochrome b562 (106 aa) extrapolated
  5ms
• B domain of protein A (58 aa) extrapolated 8ms

27
Trp-cage Folding

• Newly designed 20-residue mini-protein called
  Trp-cage
• Residue Trp6 is caged
• Folds in about 4ms
• Well defined two-state folder
• Up to 80-90 population at low temperature

Trp-cage
28
Trp-cage Folding
J. Neidigh, R. Fesinmeyer, and N. Andersen,
Nature Struc. Bio. 9, 425, 2002
29
Trp-cage Folding

• Folds in 4ms at 300K
• Well defined two state folder
• Very high population at low T

L. Qiu, S. Pabit, A. Roitberg, and S. Hagen, JACS
124, 12952, 2002
30
Trp-cage Folding
Blue MD simulation Grey NMR structure 0.97 A Ca
RMSD 1.4 A RMSD heavy atoms

• AMBER99 Force Field
• Continuum Solvent GBSA
• NVT ensemble

C. Simmerling, B. Strockbine, A. Roitberg, JACS
124, 11258, 2002
31
Trp-cage Folding

• OPLS united atom Force Field
• Continuum Solvent GBSA
• Langevin dynamics
• Water viscosity g91/ps

B MD simulation C NMR structure 2.1 A Ca
RMSD Folding time 1.5ms (3.0 A cutoff) to 8.7 ms
(2.5 A cutoff)
M. Snow, B. Zagrovic, V. Pande, JACS 124, 14548,
2002
32
Protein (un)Folding Example a b-hairpin
Protein G (2gb1)

• Res. 41-56

GEWTYDDATKTFTVTE
V. Munoz, P. Thompson, J. Hofrichter, W. Eaton,
Nature, 390, 196, 1997 R. Zhou, B. Berne and R.
Germain, PNAS, 98, 14931, 2001
33
Exploring b-hairpin Free Energy Landscape

• Implicit Solvent
• Surface Generalized-Born model
• Replica Exchange Method with 18 replicas covering
  270K to 695K
• Hybrid Monte Carlo (HMC)
• OPLSAA Force Field with SGB solvent model
• 3ns MD for each replica
• Configurations saved every 80fs

• Explicit Solvent
• 38 A cubic SPC water box (4500 atoms)
• Replica Exchange Method with 64 replicas covering
  270K to 695K
• P3ME/RESPA algorithm
• OPLSAA Force Field with SPC water
• 3ns MD for each replica
• Configurations saved every 80fs

34
(No Transcript)
35
b-hairpin Folding in Various Models

• OPLSAA/SPC (explicit)
• OPLSAA/SGB
• OPLSAA/PB

R. Zhou, B. J. Berne, PNAS 99, 2002 R. Zhou, G.
Krilov, B. J. Berne, JPC, 2004
R. Zhou, et al, PNAS 98, 2001 R. Zhou, and B.
Berne, PNAS 99, 2002
36
Lowest free energy structures
Explicit
SGB
PB

• Erroneous salt-bridges exist in all continuum
  solvent models
• Overly strongly salt-bridge effects expelled F50
  out of the hydrophobic core in SGB
• PB models behaves significantly better than the
  GB model
• Both PB and GB models need improvements

37
IV. Advanced Techniques
38
Bottlenecks in Molecular Dynamics

• Long-range electrostatic interactions O(N2)
• Intrinsic different timescales, very small time
  step needed
• Barrier crossing sampling (rare event) standard
  MD or MC not very efficient for sampling
  conformational space

39
Challenges Efforts

• Long-range electrostatics
• Particle-Particle Particle-Mesh (P3M, Hockney
  Eastwood, 1988), O(NlogN)
• Fast Multipole Method (FMM, Greengard Rhoklin,
  1990), O(N)
• Particle Mesh Ewald (PME, Dardon, 1993), O(NlogN)
• Intrinsic different time scales
• SHAKE/RATTLE (Rykaert et al. 1977, Andersen 1983)
• Constant long-range force approximation (Teleman,
  1986)
• Taylor series expansion for long-range force
  (Swindoll, 1984-86)
• RESPA (Tuckerman Berne, 1992)
• Mollified RESPA or Verlet-I (Skeel, 1997-2001)
• Barrier crossing rare event sampling
• Hybrid Monte Carlo (Duane, 1987)
• Potential smoothing with diffusion equation
  (Scheraga, 1990)
• Jump Walking (Doll Freeman, 1992-95)
• Multicanonical ensemble (Berg, 1992)
• Smart Walking (Zhou Berne, 1997)
• Replica Exchange Method (REM, Nemotto et al.
  Parisi et al., 1996-1998)
• Simulated tempering (Lyubartsev et al. Parisi et
  al., 1992)
• Catalytic tempering (Stolovitzky Berne, 2000)

40
Fast Multipole Method

• Represent charge distributions in a
  hierarchically structured multipole expansion
• Translate distant multipoles into local electric
  field
• Particles interact with local fields to count for
  the interactions from distant charges
• Short-range interactions are evaluated pairwise
  directly

CPU O(N)
41
Fast Multipole Method Procedure

• Divide the system into a hierarchical box tree
• Calculate the multipoles at lowest level directly
• Shift the multipoles from lower level to upper
  level (upward pass)
• Calculate local fields from near-distant
  multipoles (distant from the child box but not
  from the parent box)
• Shift local fields from upper level to lower
  level (downward pass)
• Evaluate potential energy and forces

42
Particle Mesh Ewald

• Similar to Ewald method except that it uses FFT
• P3ME method has a very similar spirit with PME

• Assigning charges onto grids
• Use Fast Fourier Transform to speed up the
  k-space evaluation
• Differentiation to determine forces on the grids
• Interpolating the forces on the grid back to
  particles
• Calculating the real-space potential as normal
  Ewald

43
Particle Mesh Ewald
44
(No Transcript)
45
(No Transcript)
46
Multiple Time Step Algorithm
Reversible Reference System Propagator Algorithm
(r-RESPA)
Tuckerman, Berne, Martyna, 1992
The Liouville Operator
r(tDt), v(tDt)
r(t), v(t)
The Liouville Propagator and the state of system
is given by
Trotter expansion of the Liouville Propagator
Reference System Propagator dt Correction
Propagator Dt n dt
47
RESPA for Biosystems
5-stage r-RESPA decomposition for biological
systems
Zhou Berne, 1997
The Liouville Operator decomposition for
biosystems
Reference Correction
time step
0.50 fs 1.0 fs 2.0 fs 4.0 fs 8.0 fs
48
Coupling of RESPA with FMM
5-stage r-RESPA decomposition for biological
systems
Zhou Berne, 1997
The Liouville Operator decomposition for
biosystems
Reference Correction
FMMs nearest neighbor boxes FMMs
2nd-nearest neighbor boxes FMMs multipole
expansion part
49
Coupling of Ewald, PME, P3ME with RESPA
Ewald Sum
real-space energy
k-space energy
erf(x) 1 erfc(x)
k-space interactions are not all long-range,
they still have significant contributions from
short-range while real-space interactions are
indeed short-range
PME and P3ME speeds up the k-space calculations
with FFT, O(NlogN)
50
(Ewald, PME, P3ME)/RESPA

• Normal RESPA split
• Real space short-range
• K-space long-range
• straightforward to implement
• But, k-space not all long-ranged
• Speedup is limited

• Better RESPA split
• Squeeze short-range contributions out of the
  k-space
• Modified k-space as the truly long-range force
• Modified real space as short-range force
• The modified real space can be further splited
  into two regions, short med

P. Procacci, T. Darden and M. Marchi, J . Phys.
Chem. , 1997 R. Zhou et al., J. Chem. Phys. 115,
2348, 2001
51
(Ewald, PME, P3ME)/RESPA
A better RESPA split
Using identity erf(x) erfc(x) 1
k-space energy
real-space energy

• Short-range contributions in k-space are fully
  squeezed
• k-space becomes truly long-ranged interactions
• real-space is still short-range (can be
  further separated into short med ranges)
• PME and P3ME can be implemented the same way

R. Zhou, E. Harder, H. Xu and B. Berne, J. Chem.
Phys. 115, 2348, 2001
52
FMM/RESPA Performance
53
P3ME/RESPA Performance
54
P3ME/RESPA Performance
55
Replica Exchange Method

• Jump Walking Method
• Smart Walking Method
• Replica Exchange Method

J. D. Doll and F. Freedman, J. Chem. Phys.
1992 K. Hukushima and K. Nemotto, J. Phys. Soc.
Japan, 1996 R. Zhou, B. Berne, J. Chem. Phys.
1997
56
Replica Exchange Method

• Two Step Algorithm
• (i). Each replica is simulated simultaneously
  and independently for a certain number of MD or
  MC steps
• (ii). Pick some pairs of the replicas and
  exchange the changes the configurations according
  to the following acceptance

where
57
Replica Exchange Method
58
(No Transcript)
59
(No Transcript)
60
(No Transcript)
61
Simulated Annealing
Cooling Schedules

• Most popular global optimization algorithm
• Start at high T, decrease T in small steps
  (cooling schedule)
• Easy to understand implement
• Drawback might be easily trapped in local minima

62
Umbrella Sampling
The probability density distribution in NVT
The average of a property X(p,q)
For those X(q), integrate out the momenta p, we
can have a simpler expression
B(q) any biasing function
B(q) is usually written in terms of an umbrella
potential
The umbrella potential is often harmonic, e.g.
63
Free Energy Perturbation
The Helmhotz free energy A
Then the free energy difference between two
macroscopic states is
Multiplying the numerator by a unity exp(-bH1)
exp(bH1) 1
To be accurate, we need (V2-V1) to be small, so
many stages or windows are needed. l is a
parameter, where l0 for state 1 and l1 for
state 2 (thermodynamic perturbation)
64
Further Readings Topics

• CPMD (Car-Parrinello Molecular Dynamics)
• Path Integral MD
• QM/MM methods
• Further Readings
• Computer Simulations of Liquids, Allen
  Tildesley
• Understanding Molecular Simulations, Frenkel
  Smit
• The Art of Molecular Dynamics Simulation,
  Rapaport

65
Free Energy Perturbation
Thermal Integration
The Helmhotz free energy A
Define a potential energy V(l), where l0 for
state 1 and l1 for state 2
The derivative of Helmhotz free energy A(l) with
respect to l
Then the total free energy difference between
state 1 and 2 can be integrated as