Loading in 2 Seconds...
Loading in 2 Seconds...
Molecular Biophysics III – dynamics Molecular Dynamics Simulations - 01/13/05 Hagai Meirovitch. References: A.R. Leach, “Molecular Modeling”, second edition, Prentice Hall pp. 353-406.
Molecular Dynamics Simulations - 01/13/05
Biological processes are commonly studied by experimental techniques (X-ray, NMR, etc.). However, to gain deeper insights in terms of atomic interactions we try to model biological macromolecules (proteins, DNA, carbohydrates, etc.) and simulate their behavior by Monte Carlo(MC) methods or molecular dynamics(MD) techniques that obey the rules of physics.
Before discussing MD let us refresh some basic notions from high school physics.
Examples of F (F vector; F=|F|) :
1)Stretching a spring by a distance x: F= -fx, Hook’s Law
f-spring constant. F negative because it is opposite to x.
2)Gravitation force:F= kMm/r2 - mandMmasses with distance r; k - constant.On earth(R,Mlarge),g=kM/R2 F=mg
3)Coulomb law:F=kq1q2/r2 , q1 and q2 charges.
Newton’s second law (F -resultant force):
a, v,andxvectors; ttime
If a constant force is applied along distanced,W=Fd (F=|F|).More general, W=!F.dx.
If massmis raised to height hnegative work is done, W = –mghand the mass gains potential energy,Ep= -W = +mgh- the ability to do
mechanical work: when m falls dawn, Epis converted into
Ek = mv2/2,wherev2/2=gh(at floor).
A spring stretched byd: Ep= -W = f!xdx = fd2/2
Two charges:Ep= kq1q2/r
In a closed system the total energy, Etot = Ep+Ekis constantbut Ep/Ekcan change; e.g., oscillation of a mass hung on a spring.
The total momentum of a system of particles is equal to the momentum of a single particle with the total mass of the system moving with the velocity of the center of mass of the system, vCM
vCM = Σmivi / Σmi
In a system with no external forces the total momentum is conserved - the center of mass moves in a straight line with constant speed.
Notice, the force exerted on a particle is the negative derivative of the potential energy with respect to the coordinates.
F = dp/dt = - dEp/dx
We treatN argon atoms of mass m enclosed in an isolated container. Each pair interacts via Lennard-Jones potential energy (Etot= const.)
The force in x1directionbetween (x1,y1,z1) and (x2,y2,z2)
initial velocity vi(0)(i=1,N); one can solve numerically Newton’s equations obtaining the positions xi(t) and velocities vi(t) at time t. This is the essence of MD.
Integration of the equations of motion by a finite differences method with the popular Verlet algorithm: Taylor expansions to 3rd order for i
Adding these equations gives [up to order (δt)4]
Independent of velocities. a(t) is calculated from F/m; m – mass.
The velocities (required for the kinetic energy and temperature) are:
i.e., correct up to (δt)2 . They can also be estimated at half step, t+(1/2)δt,
The process is iterative. Starts with initial coordinates and velocities; thent+δt t and t t- δt etc. Irrespective of initial conditions the particles get mixed according the laws of statistical mechanics.
Most of computer time is spent on calculation of the forces, f =a/m.
The Verlet algorithm satisfies time reversal r(t + δt) = r(t - δt).
The procedure is approximate if δt is not small enough numerical instabilities and drift in the total energy can occur. For proteins, typically δt=1/2 - 4 femtoseconds (1 fs = 10-15 seconds) large protein systems(~104atoms)are limited to ~1 ms simulations.
Many other algorithms have been developed. Some are equivalent to Verlet’s method - leap-frog (Hockney, 1970) and velocity Verlet (Swope, Andersen, Berens& Wilson, 1982), where v is more accurate
Here, first r(t+ δt) is calculated from v(t) and a(t).v(t+ δt ) is calculated in two stages, first at mid-step, i.e., v(t+ δt/2)
Finally, a(t+ δt) is calculated and the corresponding force, which lead to v(t+ δt)
The argon atoms in the isolated box (i.e., of constant energy, Etot) are described in statistical mechanics by a microcanonical ensemble; each system configuration (xN,vN ) =(x1,..,xN,v1,…,vN)has the same probability and total energy, Etot. One can calculate ensemble averages of potential & kinetic energy < Ep> + <Ek>=Etot in phase space, ΩE(tot)
One can calculate in this ensemble the distribution of the velocity component, vxof atom i (with mass m); it is the Maxwell-Boltzmann (Gaussian) probability (T – temperature; kB - Boltzmann const.)
The microcanonical ensemble provides a staticprobabilistic picture, while MD defines a deterministic dynamical picture. In other words, with MD (like in real experiments) measurements are carried out in time and it would be beneficial if the two pictures could be reconciled, i.e., if the MD time averages would lead to the ensemble averages.
Under the ergodic hypothesis a long MD run does not depend on the initial conditions and it leads to the ensemble averages. For Ep,
Where the bar denotes time average and the brackets ensemble average.
In practice, an MD trajectory after equilibration visits only a very limited part of phase space, consisting, however, of typicalcoordinates and velocities. The velocities are distributed Maxwell Boltzmann, which constitutes a criterion to check that equilibration has occurred.
A starting non-typical a typical random
configuration configuration at highT
While the probabilityof (x1,..,xN,v1,…,vN) of picture 1 is equal to that of picture 2, the number of randomly distributed configurations (as in 2) is much larger than the cluster-like conf. of fig. 1; therefore, after equilibration the system will “always” be found in a typical (random) configuration (2) and the ensemble averages will be obtained.
Calculation of temperature
In a microcanonical system(N,V,Econstant - Vvolume)T can be defined after MD equilibration, using the relation,T =<mv2>/kB>for one degree of freedomor from the total kinetic energy.
Where Nc is the number of constraints. For example, if the system is defined with periodic boundary conditions, it might be beneficial to choose initial velocities for which the velocity of the center of mass is zero. In this case Nc=3.
Periodic boundary conditions: the atom that should have left the box to position (xo+Δx,y), appears instead on the other side at position (Δx,y) with the same velocity.
The total kinetic energy is estimatedfrom the values Ek(t) obtained from a sample ofn snapshots taken at constant time intervals t:
216 argon atoms run by velocity Verlet algorithm at constant energy (NVE).
Temperature 200 K (Tcritical 150 K) – dense gas.
Time step δt= 5 femtoseconds
Each frame is taken after 20 time steps (0.1 ps)
Total run time for 2500 frames is 250 ps.
To demonstrate the independence of the equilibrium behavior on the initial conditions (ergodic theorem) the simulation starts from a non- typical configuration where the atoms are arranged at the corner with a density of liquid argon (i.e., 8 times higher than the system density after equilibration).
As we have seen, in the microcanonical ensmble (N,V,E), an MD run is completely mechanistic. However, in the canonical ensemble, where the temperature, T replaces the energy, E(i.e.,the variables are N,V,T) one has to provide a “thermostat”, i.e., a procedure to keep Tconstant. The most common thermostats are those of Berendsen and Andersen.
In the canonical ensemble the system is in contact with a large (infinite) heat bath with constant Tbath, which exchanges energy with the system. At equilibrium the average temperature of the system is also Tbath, but this temperature slightly fluctuates around its average value (Tbath), the larger the system the smaller the fluctuations. According to the Berendsen thermostat at each time step the velocities are rescaled by the factor
δt -time step
Forδt < tTthe change in velocities is more moderate and tT= 0.4 ps has been found to be appropriate.
Evaluation: easy to implement, works well, but doest not obey the canonical distribution (the velocities are not distributed according to Maxwell-Boltzmann). This procedure can lead to local hot and cold regions in the system while Tsystem is ok, i.e. Tsystem~Tbath.
With the Andersen thermostat every predefined time interval an atom is selected at random and its velocity is redefined by drawing a new velocity from a Maxwell-Boltzmann distribution. Thus, during each interval the system moves at constant energy, which is changed from interval to interval. The interval is proportional to,
is the thermal conductivity
Evaluation: The velocities are distributed as needed - according to Maxwell-Boltzmann. The trajectory is not smooth - a collection of short microcanonical stretches. If the interval is large the distribution is not canonical (close to microcanonical). If the interval is small the velocities are changed too frequently and the fluctuation of the kinetic energy is incorrect.
Methods exist where several or all the particles are treated at once.
MD not only provides statistical mechanics averages, but also allows calculating dynamical properties. One can calculate from an MD trajectory time correlation functions that lead to dynamic parameters such as diffusion coefficients. For example, the autocorrelation function of the velocity
can be estimated from an MD trajectory, where the velocities are measured n times in time intervals t
In practice, the measurements start at time (m+1)t and go back m time intervals. Contribution from all particles are considered and averaged.
In the same way one can estimate< |ri(t) - ri(0)| >2by averaging over all particles. It can be shown(Einstein, Green - Kubo) that
Where v and r denote particle velocity and position vectors. Dis the self diffusion coefficient. Notice that t should be large and in practice the length of the trajectory might be insufficient.
Other correlation functions lead to corresponding transport properties.
The typical potential energy function (force field) is:
EFF= bondsKr(r - req)2 + anglesK( - eq)2
+ dihedralsVn /2 [1 + cos(n - )]
+ i<j [Aij /Rij12 – Bij /Rij6 + qiqj /Rij]
whereKr ,K ,req,eq,,n,Vn,Aij,Bij,andqiare parameters optimized by applying this function to a large amount of experimental data and results obtained from quantum mechanical ab initio calculations.
Rijis the distance between atoms iandj and is a dielectric constant.
A protein in vacuum– no solvent effects - the screening of the Coulomb potentials by water can partially be obtained by increasing .
Typically the protein is immersed in a ‘box’ of water molecules where the system consists of n ~10,000 atoms (or more). Most of computer time is spent on calculating the forces related to water-water protein-protein and water-protein interactions.
Notice that the “spring constants” Kr andK are strong (leading to fast motions) and would require small time step, δt , relatively short trajectories. Therefore, in most studies part or all of the bond lengths are constrained by procedures such as SHAKE and RATTLE.
For proteins MD is significantly more efficient than Monte Carlo (MC) techniques, because with MD the atoms move according the calculated forces while with MC the atomic positions are changed at random which leads to high energy structures that are rejected in the process.
In an MD simulation of a protein at room temperature the system can get trapped for a long simulation time in the potential energy well of the starting structure and the move to more stable parts of conformational space is extremely slow, unreachable within a feasible simulation time.
most stable well
This problem can somewhat be resolved by methods such as parallel tempering, which is a hybrid of MD and MC.