Molecular dynamics
1 / 71

Molecular dynamics - PowerPoint PPT Presentation

  • Updated On :

Molecular dynamics ECE 697S: Topics in Computational Biology Naomi Fox March 13, 2006 Outline Why simulate motion of molecules? Energy model of conformation Two main approaches: Monte Carlo - stochastic Molecular dynamics - deterministic Applications of MD Speeding up MD

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
Download Presentation

PowerPoint Slideshow about 'Molecular dynamics' - omer

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
Molecular dynamics l.jpg

Molecular dynamics

ECE 697S: Topics in Computational Biology

Naomi Fox

March 13, 2006

Outline l.jpg

  • Why simulate motion of molecules?

  • Energy model of conformation

  • Two main approaches:

    • Monte Carlo - stochastic

    • Molecular dynamics - deterministic

  • Applications of MD

  • Speeding up MD

Why simulate motion l.jpg
Why simulate motion?

  • Predict structure

  • Understand interactions

  • Understand properties

  • Learn about normal modes of vibration

  • Design of bio-nano materials

  • Experiment on what cannot be studied experimentally

  • Obtain a movie of the interacting molecules

Structure implies function l.jpg
Structure implies function

  • central dogma: sequence to structure

  • After protein is generated, it folds without the help of external agents (except solvent)

  • Difficult problem for humans to predict 3-d structure from 1-d protein sequence

  • Organic chemists know how structure implies function for small molecular systems. Proteins are much larger, and more difficult. Open problem.

central dogma

Aquaporin simulation l.jpg
Aquaporin simulation

2004 Winner of Visualization Challenge in Science and Engineering,

Organized by the National Science Foundation and Science Magazine.

Aquaporin simulation6 l.jpg
Aquaporin simulation

Structure, Dynamics, and Function of Aquaporins

Method of molecular dynamics l.jpg
Method of Molecular Dynamics

  • Use physics to find the potential energy between all pairs of atoms.

  • Move atoms to the next state.

  • Repeat.

Energy function l.jpg
Energy Function

  • Target function that MD tries to optimize

  • Describes the interaction energies of all atoms and molecules in the system

  • Always an approximation

    • Closer to real physics --> more realistic, more computation time (I.e. smaller time steps and more interactions increase accuracy)

Scale in simulations l.jpg
Scale in Simulations

Hy = Ey



Monte Carlo

Time Scale

10-6 S




10-8 S






10-12 S

F = MA

10-10 M

10-8 M

10-6 M

10-4 M

Length Scale

Taken from Grant D. Smith

Department of Materials Science and Engineering

Department of Chemical and Fuels Engineering

University of Utah

The energy model l.jpg
The energy model

  • Proposed by Linus Pauling in the 1930s

  • Bond angles and lengths are almost always the same

  • Energy model broken up into two parts:

    • Covalent terms

      • Bond distances (1-2 interactions)

      • Bond angles (1-3)

      • Dihedral angles (1-4)

    • Non-covalent terms

      • Forces at a distance between all non-bonded atoms

The NIH Guide to Molecular Modeling

The energy equation l.jpg
The energy equation

Energy =

Stretching Energy +

Bending Energy +

Torsion Energy +

Non-Bonded Interaction Energy

These equations together with the data (parameters) required to describe the behavior of different kinds of atoms and bonds, is called a force-field.

Bond stretching energy l.jpg
Bond Stretching Energy

kb is the spring constant of the bond.

r0 is the bond length at equilibrium.

Unique kb and r0 assigned for each bond pair, i.e. C-C, O-H

Bending energy l.jpg
Bending Energy

k is the spring constant of the bend.

0 is the bond length at equilibrium.

Unique parameters for angle bending are assigned to each bonded triplet of atoms based on their types (e.g. C-C-C, C-O-C, C-C-H, etc.)

The hookeian potential l.jpg
The “Hookeian” potential

kb and k broaden or steepen the slope of the parabola.

The larger the value of k, the more energy is required to deform an angle (or bond) from its equilibrium value.

Torsion energy l.jpg
Torsion Energy

The parameters are determined from curve fitting.

Unique parameters for torsional rotation are assigned to each bonded quartet of atoms based on their types (e.g. C-C-C-C, C-O-C-N, H-C-C-H, etc.)

  • A controls the amplitude of the curve

  • n controls its periodicity

  • shifts the entire curve along the rotation angle axis ().

Torsion energy parameters l.jpg
Torsion Energy parameters

A is the amplitude.

n reflects the type symmetry in the dihedral angle.

 used to synchronize the torsional potential to the initial rotameric state of the molecule

Non bonded energy l.jpg
Non-bonded Energy

A determines the degree the attractiveness

B determines the degree of repulsion

q is the charge

A determines the degree the attractiveness

B determines the degree of repulsion

q is the charge

Energy model in action l.jpg

Stochastic method

aka Monte Carlo

Deterministic method

aka Molecular Dynamics (some ambiguity of how term is used)

Energy Model in Action

Goal: energy minimization

Simulating in a solvent l.jpg
Simulating In A Solvent

  • The smaller the system, the more particles on the surface

    • 1000 atom cubic crystal, 49% on surface

    • 10^6 atom cubic crystal, 6% on surface

  • Would like to simulate infinite bulk surrounding N-particle system

  • Two approaches:

    • Implicitly

    • Explicitly

      • Ex. Periodic boundary


Schematic representation of periodic boundary conditions.

Monte carlo l.jpg
Monte Carlo

Explore the energy surface by randomly probing the configuration space

Metropolis method (avoids local minima):

  • Specify the initial atom coordinates.

  • Select atom i randomly and move it by random displacement.

  • Calculate the change of potential energy, E corresponding to this displacement.

  • If E < 0, accept the new coordinates and go to step 2.

  • Otherwise, if E  0, select a random R in the range [0,1] and:

    • If e-E/kT < R accept and go to step 2

    • If e-E/kT R reject and go to step 2

Properties of mc simulation l.jpg
Properties of MC simulation

  • System behaves as a Markov process

    • Current state depends only on previous.

    • Converges quickly.

  • System is ergodic (any state can be reached from any other state).

  • Accepted more by physicists than by chemists.

    • Not deterministic and does not offer time evolution of the system.

Deterministic approach l.jpg
Deterministic Approach

  • Provides us with a trajectory of the system.

  • Typical simulations of small proteins including surrounding solvent in the pico-seconds.

Deterministic md methodology l.jpg
Deterministic / MD methodology

  • From atom positions, velocities, and accelerations, calculate atom positions and velocities at the next time step.

  • Integrating these infinitesimal steps yields the trajectory of the system for any desired time range.

  • There are efficient methods for integrating these elementary steps with Verlet and leapfrog algorithms being the most commonly used.

Md algorithm l.jpg
MD algorithm

{r(t+Dt), v(t+Dt)}

{r(t), v(t)}

  • Initialize system

    • Ensure particles do not overlap in initial positions (can use lattice)

    • Randomly assign velocities.

  • Move and integrate.

Leapfrog algorithm

Paper the role of molecular modeling in bionanotechnology l.jpg
Paper: The role of molecular modeling in bionanotechnology

  • Molecular modeling for r+d of bionanotechnology

  • Three case studies:

    • Polarizability of carbon nanotubes

    • Nanopores for DNA sequencing

    • Nanodiscs for experimenting on membrane proteins

Case study 1 towards an efficient description of carbon nanotube biomolecule interaction l.jpg
Case study 1:Towards an efficient description of carbon nanotube-biomolecule interaction

  • Single-walled carbon nanotubes (SWNT) can be used for in vivo biosensor applications.

  • SWNTs are highly-polarizable

  • Used MD to develop accurate description of the influence of biomolecules on nanotube electronic properties.

  • Furnishes ability to test nanotube electronic properties in realistic biological environments.

Polarizable swnt model l.jpg
Polarizable SWNT model

  • Studied water-nanotube and ion-nanotube interactions.

  • Engineers seek to control the transport of ions or charged molecules through SWNTs by means of external electric fields.

(a) Top view (b) side view of an ideal 12-section SWNT. © Potential energy profiles for same SWNT interacting with a water molecule of fixed orientation at various positions along tube axis.

Case study 2 nanopore device for high throughput sequencing of dna l.jpg
Case study 2:Nanopore device for high-throughput sequencing of DNA

  • With current technology, individual genomes determined to 99.99% accuracy within 2 months for 10 million dollars.

  • Describes a device that may provide alternative

Nanopore device l.jpg
Nanopore device

  • 2-nm-diameter pore in a 2-5nm thick synthetic membrane.

  • Sequence of a DNA molecule can be read by such a device through a semiconductor detector in the pore

    • A,C,G, and T nucleotides differ by only a few atoms

    • Essential to characterize conformation of DNA inside the pore in atomic detail

Typical translocation event l.jpg
Typical translocation event

  • A 1.4V bias applied to membrane.

  • 20 base-pair fragment of double stranded DNA placed in front of a nanopore.

  • End of DNA nearest to the pore is pulled into the pore by its charged backbone (a,b)

  • System reaches a meta-stable state (c) and translocation halts.

  • Base pairs start to split. Some freed nucleotides adhere to pore surface.

  • Voltage increased momentarily to drive system out of metastable state.

  • DNA exits pore. One of the bases holds firmly to the pore surface.

  • After 50ns, most of DNA has left pore. Nine of twenty base pairs are split.

Case study 3 nanodiscs discoidal protein lipid particles l.jpg
Case study 3:Nanodiscs: discoidal protein-lipid particles

  • Nanodiscs provide a platform in which to embed, solubilize, and study single membrane proteins.

  • Membrane proteins are difficult to study experimentally.

  • Often studied under non-native conditions, with artifacts unknown.

Md simulation of assembly of nanodiscs l.jpg
MD simulation of assembly of nanodiscs

  • Using an optimized starting ratio of protein to lipid, homogeneous nanodiscs form with discrete size and composition.

  • Predict structure of nandodiscs by simulating the experimental self-assembly.

  • To account for full process, used Marrink et-al coarse-grained lipid-water model.

    • Protein amino acid residues are represented by 2 beads

    • DPPC lipids represented by 12 beads.

    • 4 water molecules are represented 1 bead.

Conclusions l.jpg

  • Protein folds from 10 to 100s of s. Longest simulations ~1 s.

  • Can we speed up MD?

  • Would like:

    • larger systems

    • longer runs

    • more trajectories

  • approaches:

    • Mix in experimental data

    • Steered MD

    • Massive parallel computation

    • Course-grained methods

Steered MD example Bacteriorhodopsin

Steered molecular dynamics, performed by applying a series of external forces to retinal, allow one to extract retinal from bR once the Schiff base bond to Lys216 is cleaved

Future of md l.jpg
Future of MD

  • Review of deterministic method for molecular dynamics

  • Two examples of approaches for making MD ‘faster’

  • Questions

Md algorithm36 l.jpg
MD algorithm

  • Initialize system

    • Ensure particles do not overlap in initial positions (can use lattice)

    • Randomly assign velocities.

  • Compute Energy

    Etotal =

    Ebond + Etorsion + Ebond-angle + Enon-bond

  • Get velocities of atoms, then move for timestep (~1fs)

{r(t+Dt), v(t+Dt)}

{r(t), v(t)}

Two example approaches to speeding up molecular simulation l.jpg
Two example approaches to speeding up molecular simulation

  • Collect large number of trajectories using traditional MD, and see if folding rates agree with lab results.

  • Exploit geometry of the system to remove number of nodes that must be considered.

Paper absolute comparison of simulated and experimental protein folding dynamics l.jpg
Paper: Absolute comparison of simulated and experimental protein-folding dynamics

  • Statistically compared computational predictions and experimentally determined mean folding times for small protein (BBA5).

  • Used massively distributed cluster to accumulate thousands of 5 to 20 ns trajectories.

  • Found computational predictions in excellent agreement with experimentally determined mean folding times and equilibrium times.

Slide39 l.jpg

Approach to simulation l.jpg
Approach to simulation donate unused CPU cycles.

  • Simulated two BBA5 mutants

    • Considered fast-folding, mean folding time 10 s

    • 23 residues

    • Has -hairpin turn and -helix motif

  • Statistically, 10 out of 10,000 molecules will be folded after 10 ns.

  • Qualitatively reasonable but arbitrary criteria used to discriminate folded structures.

    • Each final conformation aligned to C positions of low energy BBA5 NMR structure and calculated r.m.s.d.

    • A conformation was folded if it had the helix, hairpin, and low r.m.s.d.C

Results l.jpg
Results donate unused CPU cycles.

  • Acquired a total of 32500 trajectories

  • Observed -hairpin in 1100 trajectories

  • Observed -helix in 21000 trajectories

  • Of 9000 double-mutant trajectories, 16 were folded after 20ns.

Results42 l.jpg
Results donate unused CPU cycles.

Limitations l.jpg
Limitations donate unused CPU cycles.

  • Criteria for discrimination of folded structures is arbitrary

  • Does not represent structure prediction, because relies on C positions of BBA5

  • Native state flexibility may be indication that force field is insufficiently accurate

  • Experimental folding rate uncertainty large (instrument latencies)

Conclusions44 l.jpg
Conclusions donate unused CPU cycles.

  • Contributes to body of knowledge about rapid protein dynamics.

  • Found following consistent between simulation and experimental evidence :

    • Helical structure in the unfolded state

    • Rate of helix formation

    • Rate of hairpin formation

    • Rate for folding

  • “The ability to investigate directly folding dynamics and heterogeneity with molecular dynamics will bring molecular dynamics to the same prominence in protein folding that it has already achieved in structural modeling.”

Paper constrained geometric simulation of diffusive motion in proteins l.jpg
Paper: donate unused CPU cycles.Constrained geometric simulation of diffusive motion in proteins

  • Old method

    • MD considers all 3N degrees of motion

    • Accurate trajectories for atoms only obtained with short time step (~femtosecond)

    • 10ns requires weeks of CPU time to compute

    • Most biologically relevant motions take place on milliseconds to seconds timescale

  • New method

    • Geometric simulation for proteins

Differs from md l.jpg
Differs from MD donate unused CPU cycles.

  • Not time-dependent, but distance dependent

    • Obtain collections of possible conformations, which define a trajectory that respects stereo-chemistry

    • Measure of the progression of motion is distance rather than time:

      • Mobility is monitored by RMSD from a reference conformer.

Two components to geometric simulation l.jpg
Two components to geometric simulation donate unused CPU cycles.

  • FIRST - Find which regions are rigid and which are flexible

  • FRODA - generate sterically accurate motions of a protein using rigidity information

First l.jpg
FIRST donate unused CPU cycles.

  • Floppy Inclusion and Rigid Substructure Topography

  • Given a structure (graph of atoms and covalent bonds), identify rigid and flexible regions

Froda l.jpg
FRODA donate unused CPU cycles.

  • Framework Rigidity Optimized Dynamic Algorithm

  • Finds continuous pathways as it traverses allowed regions within conformational phase space.

  • Evaluates non-covalent interaction constraints (i.e. steric clashes and hydrophobic tethers) simultaneously with the geometric constraints.

Steps l.jpg
Steps donate unused CPU cycles.

  • Perform static rigidity analysis on the protein (FIRST)

  • Perform geometric simulation (FRODA)

Rigidity analysis l.jpg
Rigidity Analysis donate unused CPU cycles.

  • Protein is a network. All covalent bond lengths and angles are constrained.

Finding rigid components l.jpg

Ring of 6 atoms and bonds donate unused CPU cycles.


Finding Rigid components

A bond-bending graph G is rigid in R3 iff

e  3n- 6for every subgraph of G having n graph vertices and e graph edges.


Bond Bending Graph:

|V|=6, |E| = 12

|E| = 3|V|-6

Constraints l.jpg
Constraints donate unused CPU cycles.

  • A subgraph is “over-constrained” if e > 3n-6

    • If you remove an edge, component is still rigid

  • A subgraph is “under-constrained” if e < 3n-6

    • Component has 3n-6 - e degrees of freedom.


e = 5 n = 4 3n - 6 = 9


e = 8 n = 4 3n - 6 = 5

Complexity of finding rigid components l.jpg
Complexity of finding rigid components donate unused CPU cycles.

  • Brute force method to find all rigid subcomponents takes exponential time (2N subgraphs in any graph).

  • FIRST uses pebble game algorithm, which takes O(N2) time, and scales linearly in practice.

    • Essence of pebble game: an algorithm for distributing the degrees of freedom belonging to the atoms (‘pebbles’) over the bonds (‘constraints’) to determine the rigidity.

    • More on pebble games at

First55 l.jpg
FIRST donate unused CPU cycles.

  • decomposes a protein down to its rigid and flexible components

Froda56 l.jpg
FRODA donate unused CPU cycles.

  • Uses geometrical constraints

  • Claim: geometric simulation is an efficient method to obtain the mobility that results from flexibility.

  • Uses creative algorithm for generating motion while enforcing constraints

    • Replace interatomic potentials by fictitious rigid bodies (the ghost templates)

    • Atoms bound to the vertices of ghost templates.

Ghost templates l.jpg
Ghost templates donate unused CPU cycles.

  • A ghost template is a rigid body that consists of set of mutually rigid atoms

  • Since they are rigid only 6 DOF

  • Atoms are bound to vertices of ghost templates

Ghost templates58 l.jpg
Ghost templates donate unused CPU cycles.

  • Defined by least-square fit of atoms

  • Vertex corresponds to specific edge

Move set l.jpg
Move set donate unused CPU cycles.

Fit templates


Fit atoms

Move set60 l.jpg
Move set donate unused CPU cycles.

The motion of an ethane molecule (a) initial atomic positions; (b) ghost templates; (c) random atomic displacement; (d) fitting of ghost templates to atoms; (e) refitting of atoms to ghost templates; (f) and (g) further iterations of (d) and (e); (h) until a new conformer is found.

Steric clashes l.jpg
Steric clashes donate unused CPU cycles.

  • Small bias to the atomic motion during the fitting (atoms to templates) process

    • Cutoff of 62% of radius

  • Sterics are always ensured (no rejections)

  • Similarly, hydrophobic tethers are ensured

    • Bias given if two atoms with hydrophobic tether move further apart than sum of radii plus .5 Å

Algorithm l.jpg
Algorithm donate unused CPU cycles.

  • Individual step ~.1 to ~.4Å RMSD

  • If algorithm does not converge in a fixed number of steps, restart from earlier conformer

Mobility results l.jpg
Mobility results donate unused CPU cycles.

50 CPU minutes!


110 residue

Above left: NMR ensemble of 20 conformers

Above right: single x-ray structure plus FRODA conformers

Right: RMSD for each residue, determined from figures above.

Directional motion l.jpg
Directional motion donate unused CPU cycles.

  • Given a starting point and ending point can we find a path?

    • Rather than randomly perturb atoms, bias towards desired structure

      • Include RMSD to target conformer in ghost/atom fitting

    • Can generate a pathway from one conformer to the other.

Directional motion results l.jpg
Directional motion results donate unused CPU cycles.


closed to occluded

Enzymatic protein

Has three very similar conformations (open, closed, and occluded), most significant differences in residues 14-24 (M-20 loop)/

Mobile loop in closed (red) and occluded (blue) positions.

Directional motion results66 l.jpg
Directional motion results donate unused CPU cycles.

From closed to occluded

Directional motion results67 l.jpg
Directional motion results donate unused CPU cycles.

Results68 l.jpg
Results donate unused CPU cycles.

  • FRODA can produce a sterically correct pathway (all motions within allowed regions).

  • Does not show time-dependency.

  • A flexible motion of a typical protein (~100 residues) investigated in 10 to 100 min.

Conclusions69 l.jpg

Massively parallel traditional MD donate unused CPU cycles.

Can find trajectories that are at most ~10 to ~20ns time duration.

Time dependent. Most accurate simulation of actual motion of each atoms

Geometric abstraction method

Can find a full sterically correct pathway from initial conformation to target conformation

Not time dependent. Not an accurate simulation of motion of each atom.


Questions l.jpg
Questions donate unused CPU cycles.

  • Can we make some relaxations to MD while minimizing loss of accuracy?

  • Can we make geometric simulation time-dependent?

  • Any other ideas for course-grained methods?

References l.jpg
References donate unused CPU cycles.

  • Understanding Molecular Simulation: from Algorithms to Applications. Daan Frenkel, Berend Smith. 2002. Academic Press. San Diego, California.

  • Deyu Lu, Aleksei Aksimentiev, AmyY. Shih, Eduardo Cruz-Chu, PeterL. Freddolino, Anton Arkhipov, and Klaus Schulten. The role of molecular modeling in bionanotechnology. Physical Biology, 3:S40-S53, 2006.

  • Snow CD, Nguyen H, Pande VS, Gruebele M. Absolute comparison of simulated and experimental protein-folding dynamics. Nature 420(6911):33-4, Nov 7 2002.

  • Stephen Wells, Scott Menor, Brandon Hespenheide, and MF Thorpe. Constrained geometric simulation of diffusive motion in proteins. Physical Biology 2:S127-S136, 2005.