## Dynamics Simulation

Molecules are dynamic, undergoing vibrations and rotations continually. The static picture of molecular structure provided by MM therefore is not realistic. Flexibility and motion are clearly important to the biological functioning of proteins and nucleic acids. These molecules are not static structures, but exhibit a variety of complex motions both in solution and in the crystalline state. The most commonly employed simulation method used to study the motion of protein and nucleic acid on the atomic level is the molecular dynamics (MolD) method (McCammon and Harvey, 1987). It is a simulation procedure consisting of the computation of the motion of atoms in a molecule according to Newton's laws of motion. The forces acting on the atoms, required to simulate their motions, are generally calculated using molecular mechanics force fields. Rather than being confined to a single low-energy conformation, MolD allows the sampling of a thermally distributed range of intramolecular conformation. Molecular dynamics calculations provide information about possible conformations, thermodynamic properties, and dynamic behavior of molecules according to Newtonian mechanics. A simulation first determines the force on each atom (F;) as the function of time, equal to the negative gradient of the potential energy (7) with respect to the position (x;) of atom i:

The acceleration at of each atom is determined by a¡ = F ¡ /m¡

The change in velocity v; is equal to the integral of acceleration over time. In MolD, one numerically and iteratively integrates the classical equations of motion for every explicit atom N in the system by marching forward in time by tiny time increments, At. A number of algorithms exist for this purpose (Brooks et al., 1988; McCammon and Harvey, 1987), and the simplest formulation is shown below:

x; (t + At) = x{ (t) + v{ (t)At v{(t + At) = v{(t) + a{(t)At = v{(t) + {F(x1 ...xN, t)/m}At The kinetic energy (K) is defined as

The total energy of the system, called the Hamiltonian (H), is the sum of the kinetic (K) and potential (V) energies:

where p is the momenta of the atoms and r is the set of Cartesian coordinates.

The time increment must be sufficiently small that errors in integrating 6N equations (3N velocities and 3N positions) are kept manageably small, as manifested by conservation of the energy. As a result, At must kept on the order of femtoseconds (10~15 s). Furthermore, because the forces F must be recalculated for every time step, MolD is a computation intensive task. Thus, the overall time scale accessible to MolD calculations is on the order of picoseconds (1ps= 1x10~i2s). The molecular simulation approximates the condition in which the total energy of the system does not change during the equilibrium simulation. One way to test for success of a dynamics simulation and the length of the time step is to determine the change in kinetic and potential energies between time steps. In the microcanonical ensemble (constant number, volume, and energy), the change in kinetic energy should be of the opposite sign and exact magnitude as the potential energy.

The MolD normally consists of three phases: heating, equilibration, and cooling. To perform MolD, the structure is submitted to a minimization procedure to relieve any strain inherent in the starting positions of the atoms. The next step is to assign velocities to all the atoms. These velocities are drawn from a low-temperature Maxwellian distribution. The system is then equilibrated by integrating the equations of motion while slowly raising the temperature and adjusting the density. The temperature is raised by increasing the velocities of all of atoms. There is a simple analytical function expressing the relationship between kinetic energies of the atoms and the temperature of the system:

T(t) = 1/{kB(3N - n)} £ m|v;| ¡=i where T(t) = temperature of the system at time t

(3N — n) = number of degrees of freedom in the system vt = velocity of atom i at time t fcB = Boltzmann constant mt = mass of atom i

N = number of atoms in the system

This process of raising the temperature of the system will cover a time interval of 10-50 ps. The period of heating to the temperature of interest is followed by a period of equilibration with no temperature changes. The stabilization period will cover another time interval of 10- 50 ps. The mean kinetic energy of the system is monitored; and when it remains constant, the system is ready for study. The structure is in an equilibrium state at the desired temperature.

The MolD experiment consists of allowing the molecular system to run free for a period of time, saving all the information about the atomic positions, velocities, and other variables as a function of time. This (voluminous) set of data is called trajectory. The length of time that can be saved during a trajectory sampling is limited by the computer time available and the speed of the computer. Once a trajectory has been calculated, all the equilibrium and dynamic properties of the system can be calculated from it. Equilibrium properties are obtained by averaging over the property during the time of the trajectory. Plots of the atomic positions as a function of time schematically depict the degree to which molecules are moving during the trajectory. The RMS fluctuations of all of the atoms in a molecule can be plotted against time to summarize the aggregate degree of fluctuation for the entire structure. The methods of MolD are becoming an important component for the study of protein structures in an effort to rationalize structural basis for protein activity and function.

Molecular dynamics simulations are efficient for searching the conformational space of medium-sized molecules, especially ligands in free and complexed states. Quenched dynamics is a combination of high-temperature molecular dynamics and energy minimization. For a conformation in a relatively deep local minimum, a room temperature molecular dynamics simulation may not overcome the barrier. To overcome barriers, conformational searches use elevated temperature (>600 K) at constant energy. To search conformational space adequately, simulations are run for 0.5-1.0 ps each at high temperature. For a better estimate of conformations the quenched dynamics should be combined with simulated annealing, which is a cooling simulation. Cooling a molecular system after heating or equilibration can (1) reduce stress on molecules caused by a simulation at elevated temperatures and take high-energy conformational states toward stable conformations and (2) overcome potential energy barriers and force a molecule into a lower energy conformation. Quenched dynamics can trap structures in local minimum. The molecular system is heated to elevated temperatures to overcome potential energy barriers and then cooled slowly to room temperature. If each structure occurs many times during the search, one is assured that the potential energy surface of that region has been adequately sought.

The molecular dynamics is useful for calculating the time-dependent properties of an isolated molecule. However, molecules in solution undergo collisions with other molecules and experience frictional forces as they move through the solvent. Langevin dynamics simulates the effect of molecular collisions and the resulting dissipation of energy that occur in real solvents without explicitly including solvent molecules by adding a frictional force (to model dissipative losses) and a random force (to model the effect of collisions) according to the Langevin equation of motion:

at = Ftm - yvt + Rt /mt where y is the friction coefficient of the solvent and R is the random force imparted to the solute atom by the solvent. The friction coefficient determines the strength of the viscous drag felt by atoms as they move through the medium, and y is the friction coefficient related to the diffusion constant (D) of the solvent by y = kBT/mD. At low values of the friction coefficient, the dynamic aspects dominate and Newtonian mechanics is recovered as y — 0. At high values of y, the random collisions dominate and the motion is diffusion-like.

Monte Carlo simulations are commonly used to compute the average ther-modynamic properties of a molecule or a molecular system especially the structure and equilibrium properties of liquids and solutions (Allen and Tildesley, 1967). They have also been used to conduct conformational searches under nonequilibrium conditions. Unlike MolD or Langevin dynamics, which calculate ensemble averages by calculating averages over time, Monte Carlo calculations evaluate ensemble averages directly by sampling configurations from the statistical ensemble. To generate trajectories that sample commonly occurring configurations, the Metropolis method (Metropolis et al., 1953) is generally employed. Thermodynamically, the probability of finding a system in a state with AE above the ground state is proportional to exp(- AE/kT). Thus, if the energy change associated with the random movement of atoms is negative, the move is accepted. If the energy change is positive, the move is accepted with probability exp(- AE/kT).