Molecular Simulations

1
Molecular Simulations
2
Molecular Simulations

• Molecular Mechanics energy minimization
• Molecular Dynamics simulation of motions
• Monte Carlo methods sampling techniques

3
Molecular Simulations

• Molecular Mechanics energy minimization
• Molecular Dynamics simulation of motions
• Monte Carlo methods sampling techniques

4
Energy Minimization
U is a function of the conformation C of the
protein. The problem of minimizing U can be
stated as finding C such that U(C) is minimum.
5
Units in force fields
Most force fields use the AKMA (Angstrom Kcal
Mol Atomic Mass Unit) unit system
Strange system!
6
Energy Minimization
Local minima A conformation X is a local
minimum if there exists a domain D in the
neighborhood of X such that for all Y?X in
D U(X) ltU(Y) Global minimum A conformation
X is a global minimum if U(X) ltU(Y) for all
conformations Y ?X
7
Notes on computing the energy
Bonded interactions are local, and therefore
their computation has a linear computational
complexity (O(N), where N is the number of atoms
in the molecule considered.
The direct computation of the non bonded
interactions involve all pairs of atoms and has
a quadratic complexity (O(N2)). This is usually
prohibitive for large molecules.
Reducing the computing time use of cutoff in UNB
8
Notes on computing the energy
Tamar Schlick, Molecular Modeling and
Simulation, Springer
9
Cutoff schemes for faster energy computation
wij weights (0lt wij lt1). Can be used to exclude
bonded terms, or to scale some interactions
(usually 1-4)
S(r) cutoff function. Three types 1)
Truncation
b
10
Cutoff schemes for faster energy computation
2. Switching
a
b
with
3. Shifting
or
b
11
Notes on computing the energy
Tamar Schlick, Molecular Modeling and
Simulation, Springer
12
The minimizers
Some definitions Gradient The gradient of a
smooth function f with continuous first and
second derivatives is defined
as Hessian The n x n symmetric matrix of
second derivatives, H(x), is called the Hessian
13
The minimizers
Minimization of a multi-variable function is
usually an iterative process, in which updates of
the state variable x are computed using the
gradient, and in some (favorable) cases the
Hessian. Iterations are stopped either when the
maximum number of steps (users input) is
reached, or when the gradient norm is below a
given threshold.
Steepest descent (SD) The simplest iteration
scheme consists of following the
steepest descent direction Usually, SD
methods leads to improvement quickly, but then
exhibit slow progress toward a solution. They
are commonly recommended for initial minimization
iterations, when the starting function and
gradient-norm values are very large.
(a sets the minimum along the line defined by
the gradient)
14
The minimizers
Conjugate gradients (CG) In each step of
conjugate gradient methods, a search vector pk
is defined by a recursive formula The
corresponding new position is found by line
minimization along pk the CG methods differ
in their definition of b - Fletcher-Reeves
- Polak-Ribiere - Hestenes-Stiefel
15
The minimizers
Newtons methods Newtons method is a popular
iterative method for finding the 0 of a
one-dimensional function
x0
x1
x2
x3
It can be adapted to the minimization of a one
dimensional function, in which case the
iteration formula is
The equivalent iterative scheme for multivariate
functions is based on
Several implementations of Newtons method exist,
that avoid computing the full Hessian matrix
quasi-Newton, truncated Newton,
adopted-basis Newton-Raphson (ABNR),
16
Molecular Simulations

• Molecular Mechanics energy minimization
• Molecular Dynamics simulation of motions
• Monte Carlo methods sampling techniques

17
What is a molecular dynamics simulation?

• Simulation that shows how the atoms in the system
  move with time
• Typically on the nanosecond timescale
• Atoms are treated like hard balls, and their
  motions are described by Newtons laws.

18
Characteristic protein motions
Periodic (harmonic)
Random (stochastic)
19
Why MD simulations?

• Link physics, chemistry and biology
• Model phenomena that cannot be observed
  experimentally
• Understand protein folding
• Access to thermodynamics quantities (free
  energies, binding energies,)

20
How do you run a MD simulation?

• Get the initial configuration
• From x-ray crystallography or NMR spectroscopy
  (PDB)
• Assign initial velocities
• At thermal equilibrium, the expected value of
  the kinetic energy of the system at temperature T
  is
• This can be obtained by assigning the velocity
  components vi from a random Gaussian distribution
• with mean 0 and standard deviation (kBT/mi)

21
How do you run a MD simulation?
For each time step Compute the force on each
atom Solve Newtons 2nd law of motion for
each atom, to get new coordinates and
velocities Store coordinates Stop
X cartesian vector of the system
M diagonal mass matrix .. means second order
differentiation with respect to time
Newtons equation cannot be solved analytically
Use stepwise numerical integration
22
What the integration algorithm does

• Advance the system by a small time step Dt during
  which forces are considered constant
• Recalculate forces and velocities
• Repeat the process
• If Dt is small enough, solution is a reasonable
  approximation

23
A widely-used algorithm Leap-frog Verlet

• Using accelerations of the current time step,
  compute the velocities at half-time step
• v(tDt/2) v(t Dt/2) a(t)Dt

v
t-Dt/2
t
tDt/2
tDt
t3Dt/2
t2Dt
24
A widely-used algorithm Leap-frog Verlet

• Using accelerations of the current time step,
  compute the velocities at half-time step
• v(tDt/2) v(t Dt/2) a(t)Dt
• Then determine positions at the next time step
• X(tDt) X(t) v(t Dt/2)Dt

v
X
t-Dt/2
t
tDt/2
tDt
t3Dt/2
t2Dt
25
A widely-used algorithm Leap-frog Verlet

• Using accelerations of the current time step,
  compute the velocities at half-time step
• v(tDt/2) v(t Dt/2) a(t)Dt
• Then determine positions at the next time step
• X(tDt) X(t) v(t Dt/2)Dt

v
v
X
t-Dt/2
t
tDt/2
tDt
t3Dt/2
t2Dt
26
Choosing a time step Dt

• Not too short so that conformations are
  efficiently sampled
• Not too long to prevent wild fluctuations or
  system blow-up
• An order of magnitude less than the fastest
  motion is ideal
• Usually bond stretching is the fastest motion
• C-H is 10fs so use time step of 1fs
• Not interested in these motions?
• Constrain these bonds and double the time step

27
Molecular Dynamics ensembles

• The method discussed above is appropriate for the
  micro-canonical ensemble constant N ( of
  particles) V (volume) and ET (total energy E
  Ekin)
• In some cases, it might be more appropriate to
  simulate under constant Temperature (T) or
  constant Pressure (P)
• Canonical ensemble NVT
• Isothermal-isobaric NPT
• Constant pressure and enthalpy NPH

How do you simulate at constant temperature and
pressure?
28
Simulating at constant T the Berendsen scheme

• Bath supplies or removes heat from the system as
  appropriate
• Exponentially scale the velocities at each time
  step by the factor ?
• where ? determines how strong the bath influences
  the system

Heat bath
T kinetic temperature
Berendsen et al. Molecular dynamics with coupling
to an external bath. J. Chem. Phys. 813684 (1984)
29
Simulating at constant P Berendsen scheme

• Couple the system to a pressure bath
• Exponentially scale the volume of the simulation
  box at each time step by a factor ?
• where ? isothermal compressibility
• ?P coupling constant

pressure bath
where
u volume xi position of particle i Fi force
on particle i
Berendsen et al. Molecular dynamics with coupling
to an external bath. J. Chem. Phys. 813684 (1984)
30
MD as a tool for minimization
Energy
Molecular dynamics uses thermal energy to explore
the energy surface
State A
State B
position
Energy minimization stops at local minima
31
Crossing energy barriers
State B
I
Energy
Position
DG
State A
A
B
time
Position
The actual transition time from A to B is very
quick (a few pico seconds). What takes time is
waiting. The average waiting time for going from
A to B can be expressed as
32
R.H. Stote et al, Biochemistry, v.43, no.24,
p.7687-7697 (2004)
33
Molecular Simulations

• Molecular Mechanics energy minimization
• Molecular Dynamics simulation of motions
• Monte Carlo methods sampling techniques

34
Monte Carlo random sampling

• A simple example
• Evaluate numerically the one-dimensional
  integral
• Instead of using classical quadrature, the
  integral can be rewritten as

ltf(x)gt denotes the unweighted average of f(x)
over a,b, and can be determined by evaluating
f(x) at a large number of x values randomly
distributed over a,b
Monte Carlo method!
35
A famous example Buffons needle problem
The probability that a needle of length L
overlaps with one of the lines, distant from each
other by D, with LD is
D
For L D
Method to estimate p numerically Throw N
needles on the floor, find needles that cross one
of the line (say C of them). An estimate of p is
Buffon, G. Editor's note concerning a lecture
given by Mr. Le Clerc de Buffon to the Royal
Academy of Sciences in Paris. Histoire de
l'Acad. Roy. des Sci., pp. 43-45, 1733. Buffon,
G. "Essai d'arithmétique morale." Histoire
naturelle, générale er particulière, Supplément
4, 46-123, 1777
36
Monte Carlo Sampling for protein structure
The probability of finding a protein
(biomolecule) with a total energy E(X) is
Partition function
Estimates of any average quantity of the form
using uniform sampling would therefore be
extremely inefficient.
Metropolis and coll. developed a method for
directly sampling according to the actual
distribution.
Metropolis et al. Equation of state calculations
by fast computing machines. J. Chem. Phys.
211087-1092 (1953)
37
Monte Carlo for the canonical ensemble

• The canonical ensemble corresponds to constant
  NVT.
• The total energy (Hamiltonian) is the sum of the
  kinetic energy and potential energy
  EEk(p)Ep(X)
• If the quantity to be measured is velocity
  independent, it is enough to consider the
  potential energy

The kinetic energy depends on the momentum p it
can be factored and canceled.
38
Monte Carlo for the canonical ensemble
Let
And let be the transition probability from
state X to state Y.
Let us suppose we carry out a large number of
Monte Carlo simulations, such that the number of
points observed in conformation X is proportional
to N(X). The transition probability must satisfy
one obvious condition it should not destroy this
equilibrium once it is reached. Metropolis
proposed to realize this using the detailed
balance condition
or
39
Monte Carlo for the canonical ensemble
There are many choices for the transition
probability that satisfy the balance condition.
The choice of Metropolis is

• The Metropolis Monte Carlo algorithm
• Select a conformation X at random. Compute its
  energy E(X)
• Generate a new trial conformation Y. Compute its
  energy E(Y)
• Accept the move from X to Y with probability
• Repeat 2 and 3.

Pick a random number RN, uniform in 0,1. If RN
lt P, accept the move.
40
Monte Carlo for the canonical ensemble

• Notes
• There are many types of Metropolis Monte Carlo
  simulations, characterized by the generation of
  the trial conformation.
• The random number generator is crucial
• Metropolis Monte Carlo simulations are used for
  finding thermodynamics quantities, for
  optimization,
• An extension of the Metropolis algorithm is often
  used for minimization the simulated annealing
  technique, where the temperature is lowered as
  the simulation evolves, in an attempt to locate
  the global minimum.