The true picture of a molecular system is far from the static, idealized image provided by molecular mechanics. The atoms and groups in molecules are in constant motion, and, especially in the case of biological macromolecules, these movements are concerted (correlated) and may be essential for biological function. The direct proof of this behavior is provided by X-ray and NMR experiments. Thermodynamic properties of molecules, especially when they form a complex or are immersed in a solvent, can not be derived from the harmonic approximation which inherently assumes small amplitude motions around a symmetric minimum. While molecular mechanics and the simulation methods described below share the same potential energy function they differ conceptually. In molecular mechanics, the information is derived from a single geometry of the molecule. Simulation methods require thousands to millions of geometries to produce meaningful results. The biological action of molecules frequently involves large amplitude motions (biochemists aften call them ``conformational transitions'' or ``changes''), resulting in abrupt changes in the geometry of the protein/ligand molecule as well as rearrangement of solvent molecules. The thermodynamic and energetic parameters of such systems can only be derived by realistic and explicit simulation, i.e., by leaving the potential energy minimum and probing the energy surface of the molecular system. Molecular simulations are a recent addition to molecular modeling systems and are still in active stages of development. They are also more complex conceptually than molecular mechanics. Recent monographs on this topic contain not only the physical and mathematical background necessary to use the method but also point to applications of the methods: Brooks et. al, (1988); McCammon & Harvey (1987), van Gunsteren & Weiner (1989).
There are essentially two approaches to performing molecular simulations: the stochastic and the deterministic. The stochastic approach, called Monte Carlo, is based on exploring the energy surface by randomly probing the geometry of the molecular system (or, in a statistical mechanics language, its configuration space). It is essentially composed of the following steps:
.
by changing the initial
coordinates at random.The most popular realization of the Monte Carlo method (though not the only one !) for molecular systems is the Metropolis method:
, 
, and 
.
corresponding
to this displacement.
accept the new coordinates and go to step 2.
,
select a random number R in the range [0,1] and:
accept the new coordinates and
go to step 2,
keep the original coordinates
and go to step 2.
As a result of a stochastic simulation, the large number of configurations (geometries) are accumulated and the potential energy function is calculated for each of them. These data are used to calculate thermodynamic properties of the system. The basic statistical mechanics formulae will be given later in this section. The Monte Carlo method is accepted more by physicists than by chemists, probably because MC is not a deterministic method and does not offer time evolution of the system in a form suitable for viewing. It does not mean however, that for deriving the thermodynamic properties of systems molecular dynamics is better. In fact many chemical problems in statistical mechanics are approached more efficiently with MC, and some (e.g., simulations of polymers chains on a lattice) can only be done efficiently with MC. Also, for Markov chains, there are efficient methods for deriving time related quantities such as relaxation times. Currently, the stronghold of MC in chemistry is in the area of simulations of liquids and solvation processes.
The deterministic approach,
called Molecular Dynamics (MD), actually simulates the time evolution of the
molecular system and provides us with the actual trajectory of the
system
. The information generated from simulation methods can in
principle be used to fully characterize the thermodynamic state of the system.
In practice, the simulations are interrupted long before there is enough
information to derive absolute values of thermodynamic functions, however
the differences between thermodynamic functions corresponding to different states
of the system are usually computed quite reliably.
In molecular dynamics, the
evolution of the molecular system
is studied as a series of snapshots taken at close time intervals
(usually of the order of femtoseconds, 1 fs = 10 
s). For large molecular
systems the computational complexity is enormous and supercomputers or
special attached processors have to be used to perform simulations spanning
long enough periods of time to be meaningful. Typical simulations of small
proteins including surrounding solvent cover the range of tens to
hundreds of picoseconds (1 ps = 10 
s), i.e., they incorporate thousands of
elementary time steps.
Based on
the potential energy function V, we can find components, 
,
of the force 
acting on an atom as:
This force results in an acceleration 
according to Newton's equation of motion:
By knowing acceleration, we can calculate the velocity of an atom in the next time step. From atom positions, velocities, and accelerations at any moment in time, we can 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. For the details of these algorithms, interested readers can consult Hermans (1985) and the references therein.
To start the MD simulation we need an initial
set of atom positions (i.e., geometry) and atom velocities.
In practice, the acceptable starting state
of the system is achieved by ``equilibration'' and ``heating'' runs prior
to the ``production'' run.
The initial positions of atoms are most often accepted from the prior geometry
optimization with molecular mechanics. Formally, such positions correspond
to the absolute zero temperature. The velocities are assigned randomly
to each atom from the Maxwell distribution for some low temperature (say
20 
K). The random assignment does not allocate correct velocities
and the system is not at thermodynamic equilibrium. To approach the equilibrium
the ``equilibration'' run is performed and the total kinetic energy (or
temperature) of the
system is monitored until it is constant. The velocities are then rescaled to
correspond to some higher temperature, i.e., the heating is performed.
Then the next equilibration run follows. The absolute temperature, T, and
atom velocities are related through the mean kinetic energy of the system:
where N denotes the number of atoms in the system, 
represents
the mass of the i-th atom,
and k is the Boltzmann constant. By multiplying all
velocities by 
we
can effectively ``heat'' the system. Heating can also be realized by
immersing the system in a ``heat bath'' which stochastically (i.e., randomly)
accelerates the atoms of the molecular system.
These cycles are repeated until the
desired temperature is achieved and at this point a ``production'' run
can commence. In the actual software, the ``heating'' and ``equilibration''
stages can be introduced in a more efficient way by assigning velocities
in such a way that ``hot spots'' (i.e., spots in which the neighboring atoms
are assigned high velocities) are avoided.
Molecular dynamics for larger molecules or systems in which solvent molecules are explicitly taken into account, is a computationally intensive task even for the most powerful supercomputers, and approximations are frequently made. The most popular is the SHAKE method which in effect freezes vibrations along covalent bonds. This method is also applied sometimes to valence angles. The major advantage of this method is not the removal of a number of degrees of freedom (i.e., independent variables) from the system, but the elimination of high frequency vibrations corresponding to ``hard'' bond stretching interactions. In simulations of biological molecules, these modes are usually of least interest, but their extraction allows us to increase the size of the time step, and in effect achieve a longer time range for simulations.
Figure 6.24: Periodic boundary conditions in molecular
dynamics.
Another approximation
is the united-atom approach where hydrogen atoms which do not
participate in hydrogen bonding are lumped
into a heavier atom to form a pseudo-atom of larger size and mass (e.g.,
a CH 
group).
Even supercomputers have their limitations and
there is always some practical limit on the
size (i.e., number of atoms) of the simulated system.
For situations involving solvent, the small volume
of the box in which the macromolecule and solvent are contained introduces
undesirable boundary effects.
In fact, the results may depend sometimes more on the size
and shape of the box than on the molecules involved.
To circumvent this limited box size difficulty, periodic
boundary conditions are used. This
idea is represented in Fig. 6.24.
In this approach, the original box containing a solute
and solvent molecules is surrounded with identical images of itself,
i.e., the positions and velocities of corresponding particles in all of the
boxes are identical.
The common approach is to use a cubic or rectangular parallelepiped box,
but other shapes are also possible (e.g., truncated octahedron).
By using this
approach, we obtain what is in effect an infinite sized system. The particle
(usually a solvent molecule) which escapes the box on the right side, enters
it on the left side, due to periodicity. Since MD simulations are usually
performed as an NVE (microcanonical) ensemble (i.e., at constant number of
particles, constant volume, and constant total energy) or an NVT (canonical)
ensemble, the volume of the
boxes does not change during simulation, and the constancy in the number
of particles is enforced by the periodicity of the lattice, e.g., a particle
leaving the box on left side, enters it on the right side. There are
also techniques for performing simulations in a NPT (isothermal-isobaric),
and NPH (isobaric-isoenthalpic) ensembles, where the pressure constancy during
simulation is achieved by squeezing or expanding box sizes. The constant
temperature is usually maintained by ``coupling the system to a heat
bath'', i.e., by adding dissipative forces (usually
Langevine, friction type forces) to the atoms of the system which as
a consequence affects their velocities.
However, each approximation has its price. In the case of periodic
boundary conditions we are actually
simulating a crystal comprised of boxes with ideally correlated atom
movements. Longer simulations will be contaminated with these artificially
correlated motions. The maximum length for the simulation, before artifacts
start to show up, can be estimated by
considering the speed of sound in water ( 
15 Å/ps at
normal conditions). This means that for a cubic cell
with a side of 60 Å, simulations
longer than 4ps will incorporate artifacts due to the presence of images.
Figure 6.25: Stochastic bondary conditions in molecular dynamics.
In the other popular approach, stochastic boundary conditions allow us to reduce the size of the system by partitioning the system into essentially two portions: a reaction zone and a reservoir region. The reaction zone is that portion of the system which we want to study and the reservoir region contains the portion which is inert and uninteresting. For example, in the case of an enzyme, the reaction zone should include the proximity of the active center, i.e., portions of protein, substrate, and molecules of solvent adjacent to the active center. The reservoir region is excluded from molecular dynamics calculations and is replaced by random forces whose mean corresponds to the temperature and pressure in the system. The reaction zone is then subdivided into a reaction region and a buffer region. The stochastic forces are only applied to atoms of the buffer region, in other words, the buffer region acts as a ``heat buffer''. There are several other approximations whose description can be found in the molecular dynamics monographs quoted here.
The straightforward result of a molecular dynamics simulation, a movie showing changing atom positions as a function of time, contains a wealth of information in itself. Viewing it may shed light on the molecular mechanisms behind the biological function of the system under study.
Molecular dynamics,
contrary to energy minimization in molecular mechanics,
is able to climb over small energy barriers and can drive the system
towards a deeper energy minimum. The use of molecular dynamics allows
probing of the potential energy surface for deeper minima in the vicinity of
the starting geometry. It
is exploited in a simulated annealing method, where the molecular system
is first heated to an artificially high temperature and then snapshots
of the trajectory at high temperature are taken as a starting state for cooling
runs. The cooled states obtained from hot frames correspond frequently to
much deeper potential energy minima than the original structure taken
for dynamics. It is particularly suitable for imposing geometrical constraints
on the system. Such constraints may be available from experimental results.
For example, the NOE constraints,
electron density maps from X-ray, etc., may be incorporated as templates
for molecular geometry. Molecular dynamics will try to
satisfy these artificially imposed constraints by jumping over shallow
potential energy minima -- behavior which is not possible in
molecular mechanics.
Since molecular dynamics is a very active field of research, new approaches are reported frequently. To facilitate the incorporation of new techniques and approaches, many molecular dynamics software packages are organized in a modular fashion. Modules communicate their data and results via disk files. In this way, the versatility and flexibility of the software is achieved, and new modules can be inserted easily between the old ones without the need of modifying existing computer code.
Statistical mechanics concepts can supply even more information from molecular dynamics runs, since they furnish prescriptions for how to calculate the macroscopic properties of the system (e.g., thermodynamic functions like free energy and entropy changes, heats of reaction, rate constant, etc.) from statistical averages of elementary steps on molecular trajectory. They also allow us to uncover a correlation between the motions of groups or fragments of the system as well as provide various distribution functions. Essentialy, all expressions derived below apply to both MD and MC, with the only difference being that the ``time step'' in MD should be changed to a ``new configuration'' for MC. MD was chosen here for illustration, since it is more often applied to chemical problems.
Three quantities form a foundation for these calculations: the partition function Z, the Boltzmann probability functional P, and the ensemble average, <A>, of a quantity A:
where 
and 
are all generalized momenta and coordinates
of the system (e.g., angular and linear velocities and internal coordinates);
is the Hamiltonian
for the system, k denotes the Boltzmann constant
( 
), and T is the temperature. The integration extends
over the whole range accessible to momenta and coordinates. For the
isochoric-isothermic system of constant NTV,
the Helmholtz free energy, A = U - TS = G - PV
(U, G, S and T are internal energy, Gibbs free energy, entropy and
absolute temperature, respectively) for the system of unit volume is given as:
The above equation simplifies significantly when cartesian coordinates
and cartesian momenta
are used, since the kinetic energy factor of the total energy represented
by the hamiltonian 
in the expression for Z reduces to
a constant depending upon the size,
temperature and volume of the system:
where the potential energy function 
depends only
on cartesian coordinates of atoms and the integral 
is called
a classical configurational partition function or a configurational integral.
Hence, the free energy, A, can in principle be calculated from the
integral of the potential energy function, 
. In practice, however,
calculation of the absolute value of the free energy is not possible, except
for very simple molecules or systems in which the motion of atoms is very
restricted (e.g., crystals). Note that cartesian coordinates span
the range 
and for systems which are diffusional in
nature (i.e., solutions, complexes, etc.) atoms can essentially occupy
any position in space and contributions from the scattered locations
to the value of the partition function are not negligible.
The situation is more optimistic with differences in the free energy
between two closely related states, since it can be assumed
that contributions from points far apart in the configurational space will
be similar, and hence will cancel each other in this case.
The free energy difference, 
,
between two states 0 and 1 is:
By multiplying the numerator of the fraction 
with the
unity written as:
the ratio becomes:
or, after including the integral from the denominator into the integral in the numerator:
i.e., the ensemble average of the exponential function involving the
difference between potential energy functions for state 0 and 1 calculated
at coordinates obtained for the state 0. To clarify the importance of this
statement, let us assume that we have two closely related states which
are described by potential functions 
and
, respectively. We perform the molecular dynamics simulation
using the function for the state 0, i.e., .
At each time step, 
, we take the current coordinates 
obtained for
system 0 and calculate the value not only of 
but
also of
,
i.e., we take atom coordinates derived for one state and calculate
the potential energy function for this state and for another state.
Then the difference
between these two values is calculated, and
the exponent of this difference taken, i.e.
. We also
calculate 
which is needed for the
Boltzmann
factor.
At the end of the calculations we combine the values from all n time steps
to calculate the ensemble average in eq. 6.50 as:
The above formula is given here only for illustration, since in practical
computations this integration is done differently for the sake of efficiency.
As was mentioned earlier, the difference between states 0 and 1 should be
very small (not more than 
kcal/mol), i.e., a perturbation.
For this reason the method is called a free energy perturbation.
If the difference between states 0 and 1 is larger, the process has to be
divided into smaller steps. In practice, the change from state 0 to 1 is
represented as a function of a suitably chosen coupling parameter
such that 
for state 0 and 
for state 1.
This parameter is incorporated into the potential energy function to allow
smooth transition between the two states. For example, if a fluorine atom
in C--F bond is being changed to a chlorine, the bond stretching term
in the potential energy function appears as:
where 
and 
are the stretching parameters
for C--F and C--Cl bonds, respectively, while 
and 
denote the optimal lengths for these bonds. Similar expressions exist
for all terms in the potential energy function and include provisions
for creating/annihilating atoms by converting them from/to dummy atoms.
Now, the path 
can be easily split into
any number of desired substates spaced close enough to satisfy the
small perturbation requirement. Assuming that the process was
divided into n equal subprocesses we have: 
,
where 
. However, in actual calculations,
the division is usually not uniform
to account for the fact that the most rapid changes in free energy occur
often at the beginning and at the end of the path. The simulation
is run for each value of providing free energy changes between
substates i and i+1 as:
i.e., the coordinates calculated at the substate 
are used
to calculate between states and
. Then the total change of free energy is calculated
as a sum of partial contributions:
The most important conclusion from the above is that the conversion between
two states does not have to be conducted along a physically meaningful
path. Obviously, atoms cannot be converted stepwise to one another, since
they are made of elementary particles. However, basic thermodynamics,
namely, the Hess's law of constant heat summation, ensures that the difference
between free energy depends only upon initial and final states.
Any continuous path
which starts at state 0 and ends at state 1,
even if it is a completely fictitious one, can be used to calculate this value.
There are other methods than perturbation of free energy to calculate the
difference in free energies between states. The most accepted one is
the thermodynamic integration method where the derivative
of free energy A versus coupling parameter is calculated:
It can be proven that:
in other words, the derivative of the free energy versus is equal to
the ensemble average of the derivative of potential energy versus ,
and hence:
In most cases, the calculations are run in both directions, i.e.,
and 
to provide the value of
hysteresis, i.e., precision of the integration process. Large
differences between and
indicate that either the simulation
was not run long enough for each of the substates, or that the number
of substates is too small. While checking the hysteresis represents
a proof of the precision of the calculations, it is not a proof of the accuracy
of the results. Errors in estimating the free energy changes may
result from an inadequate energy function and in most cases from
inadequate probing of the energy surface of the system.
Molecular dynamics, similar to molecular mechanics, can also suffer
from the ``local minimum syndrome''. If atom positions corresponding
to the minimum of energy in state 0 are substantially different from
those for state 1, the practical limitations of the length of the computer
simulation may not allow the system to explore the configuration
space corresponding to the true minimum for state 1. For this reason
the set of independent simulations, each starting at different atom
positions should be performed to assess the adequacy in sampling of the
configuration space for both states.
The methods above allow one to calculate in principle the free energies of
solvation, binding, etc. In the case of free energies of solvation, for
example, one solvent molecule is converted to a solute
molecule in a stepwise fashion,
by creating and annihilating atoms or changing their type
as illustrated by eq. 6.52, which in effect is creating
the solute molecule de novo. Binding energies could in principle be
calculated in a similar way by ``creating'' the ligand molecule
in the enzyme cavity. In practice, these calculations are very difficult,
since they involve a massive displacement of solvent, i.e., the initial
and final state differ dramatically. However, in drug design we are in most
cases
interested not in the actual free energies of solvation or binding, but
in differences between these parameters for different molecules, i.e.,
in values of 
(or 
if
an isothermic-isobaric system is being considered). The quantitative
method of
finding these parameters is often called a thermodynamic
cycle approach and is becoming a routine procedure in finding the
differences in free energy of binding for two different ligands or
the influence of mutation in the macromolecular receptor or enzyme
on binding. As an illustration, let us assume that we want to
find the difference in the free energy of binding, ( 
),
between two molecules 
and 
to a protein P.
Note that is simply related to the ratio of
equilibrium binding constants:
The thermodynamic cycle for this processes is represented in Fig. 6.26.
Figure 6.26: The thermodynamic cycle for two molecules and binding to
the same protein P and corresponding free energy changes 
.
The difference between free energies of binding for molecule and
is 
. Computation of
binding energies 
and 
is extremely difficult
since it involves bringing a molecule from a bulk solution to the
active site on the protein. This requires the removal of water
molecules from the active site and the need for very long
simulation times to reequilibrate the solvent and the protein.
Though such calculations
are possible in principle, they are still very prone to error and are
tremendously expensive. However, from Hess's law of constant heat
summation, the sum of all changes in free energy should be zero in
a closed thermodynamic cycle, i.e.,
. Hence:
Since eq. 6.50 requires that the potential energy function for
state 1 is calculated at the coordinates of state 0, the number of atoms
and bonds in both molecules must be equal.
The processes 
and 
involve annihilation/formation of atoms and cannot occur in biological systems,
however, they are much easier for calculations.
These pseudo reactions
are realized computationally by changing the type of atoms and groups.
If molecules are of similar size the solvent molecules need not
diffuse out of the active site. If new
atoms are needed, dummy atoms are converted stepwise to real atoms.
On the other hand, if molecule has fewer atoms than molecule ,
atoms are annihilated by converting them to dummy atoms.
As an example consider two molecules which
differ by the presence of a methyl group. ``Mutating'' the hydrogen atom
into a methyl group is illustrated in Fig. 6.27.
Figure 6.27: Conversion of a hydrogen atom into a methyl
group. Note that
the number of atoms and bonds is conserved during the conversion.
A molecular dynamics simulation is first performed to calculate
, i.e., the free energy change corresponding to
``mutating'' molecule
into in the solution. In this simulation the protein molecule
is not present in the system since it is inert in this reaction
(i.e., ``reduces'' on both sides of the chemical equation). The second
simulation involves calculation of 
, i.e., the molecule
is converted to molecule inside the active site of the protein P.
The 
is then calculated as 
.
It cannot be stressed strongly enough that both of these free energy changes,
and , are
non-physical. is not a difference between the free energy of
solvation of molecules and unless the ``molecules'' are
mononuclear. Calculation of the difference in solvation energies would
require a dynamics simulation for molecules in a vacuum to
account for the differences in the interaction energy and the entropy of
free and solvated molecules. Similarly, alone is not a difference
between the free energies of binding,
since it contains no information about that
solvent interaction with the part of the molecule which is
submerged in the active site.
The differences in the solvent-solute interaction
between ligands may comprise a major portion of their relative affinity for
the protein.
In these instances where we know the structure of the receptor protein, the thermodynamic cycle method is a very promising tool for drug design, since it allows us to actually calculate the differences between free energies of binding for molecules which have not even been synthetized yet! Likewise, it allows assessing the impact of variations in protein structure (e.g., mutant versus wild type protein) on the binding, i.e., perform genetic engineering on the computer. This is why the X-ray crystallography or NMR studies of proteins combined with molecular simulation techniques form a framework for the drug design of the future, where prospective drug molecules will be first inferred from theory, and only later synthetized. Morerover, intensive research on deducing protein tertiary structure from the knowledge of its primary sequence is bound to result in new possibilities for folding of proteins in the computer. With the astonishing advances in the power of computers and progress in methodology, there is no doubt that the importance of these techniques to rational drug design will increase steadily.
