0
Methods for
the Generation and Analysis of
Protein Simulations

Mark Gerstein
<mbg@hyper.stanford.edu>
Stanford Structural Biology 

http://hyper.stanford.edu/~mbg/SurfVolTalk/

N.B.  What follows are transparencies created for a talk. Feel free to download and read them as you would the Xerox of a paper. However, some of transparencies contain copyrighted illustrations, and others are property of the author. Consequently, you should  NOT incorporate any of them into published documents or distribute them without contacting the author.   

1
Methods for
the Generation and Analysis of
Protein Simulations

Mark Gerstein

http://hyper.stanford.edu/~mbg/SurfVolTalk/

2
Overview
1	Simulation Methods
Potential Functions
Minimization
Molecular Dynamics
Monte Carlo
Simulated Annealing 
2	Types of Analysis
liquids: RDFs, Diffusion constants
proteins: RMS, Volumes, Surfaces 

•	Established Techniques
(chemisty, biology, physics)
•	Focus on simple systems first (liquids). Then explain how extended to proteins. 
3
Potential Functions
Each atom is a point mass (m and x)
Atoms interact through a variety of forces
Also, for proteins there some special pseudo-forces: torsions and improper torsions, H-bonds. 
4
Minimization
•	Particles on an “energy landscape.” Search for minimum energy configuration
•	Steepest descent minimization
°	Follow gradient of energy straight downhill.
°	i.e. Follow the force: step ~ F = -— U
so     
•	Other methods
°	conjugate gradient
step ~ F(t) - bF(t-1)
°	Newton-Raphson: using 2nd derivative, find minimum assuming it is parabolic
•	Get stuck in local minima

5
Molecular Dynamics
•	Give each atoms a velocity.
°	If no forces, new position of atom would be determined only by velocity
x(t+∆t) = x(t) + v∆t

•	Forces change the velocity, complicating things immensely

6
Molecular Dynamics (cont)
•	On computer make very small steps so force is nearly constant and velocity change can be calculated (uniform a)
°	 
°	[Average velocity over ∆t] = (v + F(∆v,2))
°	Trivial to update positions:

•	Step must be very small
°	∆t ~ 1fs 
(atom moves 1/500 of its diameter)
°	This is why you need fast computers
•	Actual integration schemes slightly more complicated
°	Verlet (explicit half-step)
°	Beeman, Gear
(higher order terms than acceleration)
7
Phase Space Walk
•	Trajectories of all the particles traverses space of all possible configuration and velocity states (phase space)
•	Ergodic Assumption: 
Eventually, trajectory visits every state in phase space
•	Boltzmann weighting: 
Throughout, trajectory samples states fairly in terms of system’s energy levels
°	More time in low-U than high-U states 
°	Probability of being in a state ~ exp(F(-U,kT))
•	Consequently, statistics (average properties) over trajectory are thermodynamically correct
8
Example Phase Space Walk

9
Monte Carlo
•	Other ways than MD to sample states fairly and compute correctly weighted averages? Yes, using Monte Carlo calculations. 
•	Basic Idea: 
Move through states randomly, accepting or rejecting them so one gets a correct “Boltzmann weighting” 
•	Formalism:
°	System described by a probability distribution r(n) for it to be in each state n
°	Random (“Markov”) process p operates on the system and changes distribution amongst states to pr(n)
°	At equilbrium original distribution and  new distribution have to be same as Boltzmann distribution 

10
Monte Carlo (cont)
•	Metropolis Rule (for specifying p )
1	Make a random move to a particle and calculate the energy change ∆U
2	∆U < 0 Æ accept the move
3	Otherwise, compute a random number R between 0 and 1: 
R < ~ exp(F(-U,kT)) Æ accept the move
otherwise Æ reject the move

•	“Fun” example of MC Integration
°	Particle in empty
box of side 2r (energy of all states same)
°	p = 6 ¥ [Fraction of times particles is within r of center]		

11
MC vs/+ MD
•	MD usually used for proteins. Difficult to make moves with complicated chain.
•	MC often used for liquids. Can be made into a very efficient sampler.
•	Hybrid approaches (Brownian dynamics)
•	Simulated Annealing. Heat simulation up to high T then gradually cool and minimize to find global minimum. 
•	Periodic Boundary Conditions

12
Typical Systems: protein + water

[prot-wat-p-.rgb]
[prot-wat-w-.rgb]
[prot-wat-pw.rgb]
13
Average over simulation 
•	Deceptive Instantaneous Snapshots   (almost anything can happen)
•	Simple thermodynamic averages 
°	Average potential energy <U>
°	T ~ <Kinetic Energy> = 
•	Some quantities fixed, some fluctuate
in different ensembles
°	NVE protein MD (“microcanonical”)
°	NVT liquid MC (“canonical”)
°	NPT more like the real world
14
Timescales 
Motion	length	time	
	(Å)	(fs)	
bond vibration	0.1	10	
water hindered rotation	0.5	1000	
surface sidechain rotation	5	105	
water diffusive motion	4	105	
buried sidechain libration	0.5	105	
hinge bending of chain	3 	106	
buried sidechain rotation	5	1013	
allosteric transition	3	1013	
local denaturation	7	1014	
(From McCammon & Harvey, Eisenberg & Kauzmann)

15
D & RMS 
•	Diffusion constant
°	Measures average rate of increase in variance of position of the particles 
°	
°	Suitable for liquids, not really for proteins
•	RMS more suitable to proteins
°	
°	di = Difference in position of protein atom at t from the initial position, after structures have been optimally rotated translated to minimize RMS(t) 
•	Solution of optimal rotation has been solved a number of ways (Kabsch, SVD)
[draw?]
16
Number Density 
=	Number of atoms per unit volume averaged over simulation divided by the number you expect to have in the same volume of an ideal “gas”
	
•	Spatially average over all directions gives 1D RDF
=	F(Average number of neighbors at r, Expected number of neighbors at r)
°	“at r” = means contained in a thin shell of thickness dr and radius r.
 
17
Number Density (cont)
•	Advantages: Intuitive, Relates to scattering expts
•	Not applicable to real proteins
°	1D RDF not structural
°	2D proj. only useful with “toy” systems 
•	Number densities measure spatial correlations, not packing
°	Low value does not imply cavities
°	Complicated by asymmetric molecules
°	How things pack and fit is property of instantaneous structure — not average

18
Packing

close-packed protein

“Open” structure of water
[Color-Image 3a: hoh-v-ar.rgb]
[Color-Image 3b: water6ps.c.rgb.gz]
[Color-Image 3c: ar2.g.rgb.gz]
19
Voronoi Volumes
•	Applicable to the simulations of real proteins
•	Each atom is surrounded by a single convex polyhedron and allocated the space within it
°	complete allocation (no cavities)

•	Ideal for studying packing
Packing efficiency = 
20
Voronoi History
•	Mathematician Voronoi 1908
•	Bernal & Finney 1967, liquid simulation
•	Common in CS
°	“Nearest-neighbor” problems
(find the neighbor of a query point)
°	Largest empty circle in a collection of points
°	Dual of a Voronoi diagram is Delaunay triangulation (Convex Hull)
21
Basic Voronoi
     
1	Find Vertices
°	Point in space assigned to closest atom
°	In 3d, points equidistant from 2 atoms form a face; 4 atoms, a vertex
°	Solve for vertex using 4 atoms and equation of sphere. (Make sure vertex is closer to these 4 than any others)
2	Traverse Vertices Æ Calculate Volume
°	Label Vertices: those with 2 atoms in common are a face; those with 3, connected by edge 
°	Systematically traverse: face areas Æ pyramid volumes Æ polyhedra volumes
22
Problem of differently sized protein atoms
 •	Construction based around partitioning space amongst “equal” points, but not all protein atoms are equal: some are larger (e.g. S v. O)
•	Richards proposed solution in 1974
(+  Finney, 1982)

Dotted Line — Basic Voronoi bisection
Solid Line — Richards’ modification
23
Protein Surfaces
•	Why Useful?
°	Protein is solid object. Surface is where action takes place.
°	Surface useful for docking and drug-design
°	Hydrophobic energy proportional to surface area
•	Voronoi surfaces
°	Area of a face “weights” extent of contact
°	Volume as weighted average of contact distances
°	Precise calculation of contacts
•	Problem with Voronoi construction
°	Missing neighbors lead to misshapen polyhedra
°	Waters often missing in crystal structures — but not in simulation
24
Accessible Surface
•	Accessible Surface
Roll sphere (water) on surface and look at locus of sphere centers. Usually represented as a dot surface. 
•	Problem with accessible surface: It has sharp cusps. Not smooth and continuously differentiable (relevant for energy calculations)
25
Accessible Surface Calculation 
•	Lee & Richards algorithm (1st, 1970)
°	Pick an arbitrary direction from which to view the protein. Slice it into many sections perpendicular to this direction. In each section, cycle over all the atoms. Each atom is represented as a sphere with a radius that is the sum of its VDW radius plus that of a probe solvent -- i.e. 1.4 for water. For each atom determine the circle corresponding to the intersection of this sphere with the sectioning plane. Remove all parts (i.e. arcs) of this circle occluded by the circles of other atoms. Multiply the total amount of non-occluded arc length by the sectioning width to get the surface area for atom. Sum over all atoms and all sections to get total area.
•	Shrake & Rupley algorithm 
°	Surround each atom with sphere of uniformly spaced dots (e.g. 92). Remove dots contained in other atoms spheres. Total number of remaining dots is accessible surface. 
26
Molecular Surface
•	Molecular Surface
M.S. = contact surface + re-entrant surface 
C.S. = points of tangency between probe sphere and protein when probe sphere is only touching one atom
R.S. = solid angle of probe sphere when tangent to two protein atoms
•	First proposed by Richards, but hard to calculate. First numeric calc. by Connelly. Later analytic calculation by Connelly. Analytic version is continuously differentiable.

[talk.ms.gif]
27
Hydration Surface
•	Problem with molecular surface:
°	Water does not really roll on protein surface like a sphere.
°	This is treating water chemically like liquid argon.
•	Solution: calculate actual locus of real water positions in a simulation and fit a surface through the mean position of the second hydration shell. 

28
References
Allen, M. P. & Tildesley, D. J. (1987). Computer Simulation of Liquids. Claredon Press, Oxford.
Branden, C. & Tooze, J. (1991). Introduction to Protein Structure. Garland Publishing Incorporated, New York.
Daggett, V. & Levitt, M. (1993). Realistic Simulations of Native-Protein Dynamics in Solution and Beyond. Ann. Rev. Biophys. Biomol. Struct. 22, 353-380. 
Eisenberg, D. & Kauzmann, W. (1969). The Structure and Properties of Water. Clarendon Press, Oxford.
Eisenberg, D. & McLachlan, A. D. (1987). Solvation energy in protein folding and binding. Nature 319, 199-203. 
Gerstein, M. & Lynden-Bell, R. M. (1993a). Simulation of Water around a Model Protein Helix. 2. The Relative Contributions of Packing, Hydrophobicity, and Hydrogen Bonding. J. Phys. Chem. 97, 2991-2999. 
Gerstein, M. & Lynden-Bell, R. M. (1993b). What is the natural boundary for a protein in solution? J. Mol. Biol. 230, 641-650. 
Gerstein, M., Tsai, J. & Levitt, M. (1995). The volume of atoms on the protein surface: Calculated from simulation, using Voronoi polyhedra. J. Mol. Biol. 249, 955-966. 
Lee, B. & Richards, F. M. (1971). The Interpretation of Protein Structures: Estimation of Static Accessibility. J. Mol. Biol. 55, 379-400.
Levitt, M. (1982). Protein conformation, dynamics, and folding by computer simulation. Ann. Rev. Biophys. Bioeng. 11, 251-271. 
Levitt, M. & Sharon, R. (1988). Accurate Simulation of Protein Dynamics in Solution. Proc. Natl. Acad. Sci. USA 85, 7557-7561. 
McCammon, J. A. & Harvey, S. C. (1987). Dynamics of Proteins and Nucleic Acids. Cambridge UP, 
Richards, F. M. (1974). The Interpretation of Protein Structures: Total Volume, Group Volume Distributions and Packing Density. J. Mol. Biol. 82, 1-14. 
Richards, F. M. (1977). Areas, Volumes, Packing, and Protein Structure. Ann. Rev. Biophys. Bioeng. 6, 151-76. 
Richards, F. M. (1985). Calculation of Molecular Volumes and Areas for Structures ofKnown Geometry. Methods in Enzymology 115, 440-464. 
Shrake, A. & Rupley, J. A. (1973). J. Mol. Biol. 79, 351. 

More at
http://hyper.stanford.edu/~mbg/SurfVolTalk/