Molecular Mechanics Minimization:
 Minimization
 Exploring Multidimensional PES

MD:
 MC
 Analysis:
Minimization
Object: To find a set of coordinates corresponding to minima on the
PES.
Requires: Taking the derivative of the potential function w.r.t.
all d.o.f. in the molecule, and locating where each of those derivatives
are simultaneously zero.
Verify:. The hessian is used to verify the minimum
conformations
Why minimize?
 Comparison of structures/properties
 Preparation for MD
Since one begins with zero force in MD, a zero E state is desired in the
initial step. Thus minimization is first performed on the molecule.
 Template forcing
 Systematic mapping of E space
 Binding energies
 Docking
 Harmonic analysis
 Comparing/Fitting force fields
 Algorithms:
There are a plethora of algorithms for performing function
minimization. Although many of these procedures are adequate for minimizing
5 10 variables most are not well suited for MM calculation which typically
involve over 100 atoms.
Methods are either 1) convergent  which generally stop when they have
achieved a minimum energy conformation or, 2) nonconvergent  which never
stop wandering on a PES.
 Steepest Descents
 Conjugatedirection
 Newton methods
 QuasiNewtonRaphson
 Full NewtonRaphson
 Quenched dynamics
Minimum information to start minimization: Starting set of
atomic coordinates
 Parameters for various terms of the potential energy function
 Description of molecular topology (covalent structure) (potential
energy function)
 Additional features
 Use of periodic boundary conditions so that the simulated system is
effectively part of an infinite system.
 Forces of constraint to hold the system in a particular configuration
or move it toward a different target conformation.
 Provision for occasionally adjusting positions, velocities or
parameters in order to maintain or change T,P,V, forces of constraint or
other externally controlled conditions.
Energy minimization Algorithms:
Energy, first and second derivative used for any set of atomic positions to
generate new sets of coordinates in effort to reduce total P.E. and, by
repeating, to optimize.
Application in Macromolecular Structure:
 Assist in refinement:
 Structures of peptides (Scheraga, 1968)
 Crystal structure of proteins (Levitt and Lifson, 1969; Levitt, 1974)
 Nucleic Acids (Jack et al., 1976)
 Used to optimize bond lengths and angles in crystallographic refinement
procedures. (Jack and Levitt, 1978; Konnert and Hendrickson, 1980)
 In addition to problems of static nature, provide information on
dynamic processes.
If a macromolecule is deformed in steps along a chosen reaction coordinate
and the rest of the structure relaxed by minimization while the reaction
coordinate is held fixed, the variation in conformational energy along
reaction coordinate gives an estimate of the shape of the energy barrier
(basis for adiabatic mapping).
 Molecular Dynamics
Minimization and Adiabatic Mapping:
 When used to examine structural transitions the energies that are
calculated are only approximations to the familiar thermodynamic functions.
In the absence of thermal motion, there is no entropic contribution to the
calculated energy. The calculated energy is more closely related to the
thermodynamic internal energy or enthalpy than to a free energy, but does
not include the T dependent contributions to these functions.
 For static models, the vibrational entropy and other thermal
contributions can be estimated from normal mode analysis.
 MD, on the other hand, explores the confrontational space in a fashion
that directly reflects the density of states, so entropic and other thermal
effects are naturally included. Therefore, free energy differences can be
determined by a molecular dynamics simulations of different energy
conformations along a transition pathway.
Minimization:
 Nonlinear optimization
 Given a set of independent variables and a specific objective function,
the task is to find a set of values for the independent variables for which
the function has its minimum value.
 In the case of a model macromolecule with N atoms, the 3N components
are he atomic coordinates and the objective function is the potential
energy function.
 Macromolecule Modeling:
 Steepest Descents
 Conjugate Gradients
 Smaller Molecules
 Steepest Descents
 Conjugate Gradients
 Newton Raphson
 Improving Computational Efficiency
I. Nonbonded Cutoffs
 Several approximations can be made to the force field to improve
the computational efficiency. The most common of these is a nonbond cutoff
(i.e., neglecting nonbonded interactions for pairs of atoms separated by
distances greater than a cutoff value.
 This can be a rather drastic simplification.
See review Brooks, et al. (1985)
Figure: Proteinsubstratesolvent system with 5000 atoms
(typical of a small protein (100150 residues) surrounded by a 12 ;layers
of water. Figure shows the number of nonbonded pairwise interactions in
millions expected for a 5000 atom system as a function of cutoff distance.
The time required to evaluate the total energy of this system is
approximately proportional to the number of nonbonded interactions.
The calculation would run at least 10 times faster with an 8 A cutoff than
with no cutoff (assuming the usual case that the nonbond term is rate
limiting).
 The tradeoff is that interactions beyond the cutoff distance are not
accounted for. Both van der Waals interactions and electrostatic
interactions are significant up to 15A.
Figure: Calculation of the nonbonded energy as a function of
cutoff distance in the [AlaProDPhe]2 crystal. Kitson and Hagler have
shown that the energy accounted for changes from 63% to 97% of the
asymptotic value as the cutoff distance is increased from 8 to 15 A (1986).
This figure shows how the van der Waals component of the nonbond energy
varies as a function of cutoff distance for this system. The exact
dependence of the energy on the cutoff distance depends on the system
itself and should be calibrated for each new system of interest.
II. Modeling and Design Strategies using Minimization:
 A. Constraints to the Target Function
 B. Restraints to the Target Function
 C. Exploring Conformational Space
A. Constraints
 Modifying the target function is a powerful way to customize a
minimization to address specific modeling objectives.

We define constraints as d.o.f. which are fixed, i.e., not allowed to vary
during the course of minimization.
 Impact on computation efficiency:
 When constraints are in place, the target function may then have
fewer terms to evaluate and can usually be evaluated faster,
 Since the number of iterations required to converge is roughly
proportional to the number of d.o.f. in most algorithms, reducing the
degrees of freedom translates directly into faster minimizations.
The tradeoff is that the resulting model may include unrealistic strain of
other artifacts that could compromise structural and energetic analysis.
Use: In preliminary studies to quickly setup systems for more
subsequent more rigorous calculations
 When computer time for unconstrained calculations is unavailable.
Types:
 Rigid body minimization
 Torsional minimization

Fixed atom minimization
1.Rigid body minimization
 If all internal motions are constrained, translations and rotations
are the only remaining d.o.f.
 Uses:
 Useful for studying the docking of two(or more) relatively
rigid molecules.
 For molecules that are not particularly rigid, such as the active site
of an enzyme rigid body dockings can quickly position the substrate in a
likely orientation in preparation for more rigorous minimizations or
dynamics simulation.
 Studying macromolecular packing.
Advantages:
 Minimizations can be performed quickly, even to the
point of their being interactive (docking substrates into proteins) (target
function has fewer terms to evaluate)

Qualitative assessments before full minimizations can be performed.
Restriction:
No way for the internal d.o.f. to absorb any strain
introduced in the docking process.
2. Torsional minimization
 Since bond length and angle distortions in general require
significant energy, most bond lengths and angles in amino acids vary only
slightly from average values. Thus, to a first approximation, the
confrontational flexibility of proteins can be attributed to low energy
rotations about single bonds. By only allowing motion about dihedral
angles, the number of d.o.f. can be reduced from 3N to approximately N/2
(the number of nonhydrogen containing bonds) or even lower if bonds with
double bond character are constrained.
 Well known quantitative problems:
Flexible geometry map has a
larger region of accessible conformational space, since it allows bad
contacts resulting from torsional displacements to be reduced by small
length and bond angle displacements, as would be in nature.
Figure:
Restricted flexibility has quantitative effects on both entropy and enthalpy.
Figure: Shows energy barriers to rotation are consistently overestimated by
414 kcal.mole by rigid torsional minimizes for a particular molecule.
NOT appropriate to use rigid body techniques to estimate energy barriers
about bonds except for the most qualitative cases. Because torsonial
minimizers cannot distribute strain into the bond angles of the system,
they tend to overestimate the height of energy barriers and distort the
energy surface.
None of the commonly used methods generate large collective displacements,
since they do not efficiently simulate the transmission of forces through
the structure. The following two methods allow such collective
displacements.
 Pseudodihedral Method: Harvey and McCammon, 1982
 Molecule is
divided into regions by connecting pairs of atoms with virtual bonds which
serve as the axes about which rigid rotations of group of atoms takes place
 By keeping the number of virtual bonds small (<12) it becomes
feasible to perform a grid search on the conformational spaces defined by
these pseudohedral rotations.
 This local grid search method provides an initial conformation for a
more serious minimization.

Example: Pseudodihedral algorithm@ steepest descents reduced the
required computation time by ~40% Ct RNA
 Variable Metric Method plus soft potentials
 Suitable for extensive conformational changes
 Developed by
Levitt(1983) for folding model proteins into trial conformations suitable
for further refinement by more conventional methods.
 Minimizer: Variable
Metric
Potential Function: Soft atom allow atoms to pass through one another
 Minimization with fixed atoms
 The simplest constraint conceptually is to fix atoms in cartesian
space. With every constrained atom, 3 d.o.f. are removed. Depending on the
modeling objective, holding a protein fixed (or at least the part of it
that is far from the region of interest.
 Hagler/Moult, 1978
 Moult/James, 1986
 Fine et al, 1986
B. Restraints
A term added to the target function that biases the system toward adopting
a certain value for a d.o.f. (but not requiring it absolutely as in the
case of a constraint).
 No Computational advantage:

The number of d.o.f. are not reduced
The target function is actually more complicated (though not significantly)
Advantages:
 More control over the minimization
 Exploring defined regions of
the conformational space
 Testing conformational hypothesis (template forcing;Brooks et al.,1983)
Types: Torsional Restraints
 Distance Restraints
 Template Forcing
 Tethering
1.Torsonial Restraints:
 Places a harmonic torque about a bond to force it to a new value
Etorque=ktorque(00target)2
 By systematically changing the target angle and minimizing the entire
structure at each step, one can produce an adiabatic surface. These maps
are useful to quantify low energy pathways across a conformational barrier.
 The torsion angles implied by NMR coupling constants can also be used
as target angles using minimization to generate a conformation biased
toward a structure consistent with experimental data.
2. Distance Restraints
 A simple but powerful restraint is to add a term restraining the
distance between two atoms. Distance restraints are useful for
incorporating NOE data from NMR experiments into a structural model.
Including NOE distance restraints with the energyunction can produce a
family of representative structures of reasonable energy consistent with
the given data and constitutes a powerful method of solving solution
structures.
 Various functional forms are used. A harmonic function amounts to
placing a classical spring between atom pairs (Braun and Go, 1985; Van
Gunsteren and Karplus, 1980).
E=k(RijRtarget)2
 Often, one sided restraints are used. By
turning off the restraint distances greater than the target distance, one
can push the atoms apart. Conversely, if the restraint is turned off at a
distance less than the target, the atoms are pulled together.
 Constant rather than harmonic forces may be used when initial
structures are far from the target value. For example, a constant force
would gently close a ring when a harmonic force could add enough strain to
invert chiral centers.
 Most applications of NOEtype distance restraints to structure solving
have used either distance geometry algorithms
Crippen, 1977; Havel et al., 1983
or molecular dynamics
Clore et al., 1986; Havel and Wuthrich, 1984
to satisfy the initial restraints. Minimization is used only in the final
refinement stages.
 Braun and Go, 1985
Other applications have used minimization exclusively in conjunction with
perturbations to the target function to solve NOE restrained structures.
The essential feature was that long range nonbond interactions were ignored
during the initial stages of restrained minimization. This allows atoms and
residues to pass through each other as they seek to satisfy the
experimental restraints. Gradually the sphere of nonbonded interactions is
increased as the model becomes more refined. In addition, the minimization
occurred in torsion space, which accelerates convergence and allows for
larger scale motions.
3. Template forcing
 When testing conformational hypotheses for 'binding conformations in
drug design, it is often useful to know if one molecule can adopt the
conformation ( or more generally, the shape) of another.
Brooks, et al., 19832
 By selecting corresponding atoms or functional groups between two
molecules atoms or functional groups between two molecules, one can force
atoms from a flexible molecule to superimpose onto atoms of a rigid
"template" molecule. Effectively, a spring is placed between the paired
sets of atoms, pulling the atoms in the flexible molecule toward a rigid
molecule.
Etemplate=Ektemplate(Ri analog  Ri template)2
The restraining atom pairs or forcing term is added to the energy to
construct the target function to be minimized. By varying the spring force
constant, one can force the flexible molecule into the conformation of
template to an arbitrarily small deviation at the lowest energy cost. The
energy required to achieve a give 'fit' can then be used to evaluate how
easily an analog can adopt the conformation of a given template.
4. Tethering atoms to points in space.
 Special case of template forcing in which the 'template' is a set of
points in space, usually corresponding to a desired target conformation
such as a crystal structure.
 Useful for 'annealing' crude complex structures before carrying out
unrestrained minimization or dynamics.
 For example, when minimizing solvent with a protein, one can tether the
oxygens of the water to their initial positions while the hydrogens
reorient (Roberts,et al., 1986). One may gradually reduce the force
constant during the minimization and update the tethering points
periodically to correspond to the current structure.
 Warme and Scheraga, 1974
First used tethering to keep lysozyme close to the crystal coordinates
during minimization. This is desirable since the crystal data may include
disordered or poorly ordered regions resulting in strained geometrics and
overlapping atoms that would result in unrealistic long range disruptions
to the structure.
Classical Minimizers:
 Locate 'local', not 'global' minimizers.
 Designed to 'ignore' configurations if E increases
 Do not push systems over barriers, but, down to valleys.
C. Searching Conformational Space
 So far we have looked at classical minimizers (SD,CG,NR). A common
limitation of classical minimization algorithms is that they locate a local
minimum usually close to the starting configuration, and not the global
minimum. This is because these minimizers are specifically designed to
ignore configuration if the energy of that configuration increases, and
thus, do not push systems over barriers, but down valleys. Modified
strategies are needed to search conformational space more thoroughly.
 Two approaches to searching conformational space:
 Screen from a large number of systematically generated points in
conformational space, for the 'best' starting points from which to begin a
minimization.
 Modify target function and/or minimization algorithms to allow systems
to overcome local barriers enroute to finding lower energy minima.
Methods:
 Systematic searching
 "Uphill minimizers
 Minimization through reduction of dimensionality
 Dynamic searching
 Monte Carlo
1. Systematic Searching of Conformational Space:
 Generated by rotating about all torsion angles in small increments
 Only method that can ensure finding global minimum
 Severe restriction:# of calculations is an exponential function of d.o.f.
 Modifications:
 Crippen & Scheraga, 1971; 1973
Employ pruning algorithms that ignore large volumes of unlikely
conformational space.
 Scheraga; Vasquez & Scheraga, 1985; Gibson & Scheraga, 1986 "Build up
method"applied to peptides. Uses standard geometrics of smaller peptide
segments to build up dipeptides and after cutting out high energy
structures, these are used to build all possible tetrapeptides, etc. in an
iterative manner.
2. "Uphill" minimizers:
Overcomes the natural tendency of minimizers to get stuck in local wells.
 Methods:

Deflation Method: 'Synthetic Division'
Minimizes a target function that is constructed by dividing the potential E
by the found "roots" of the E surface each time a local minimum is found
and uses the new target function o continue the search.
 Gentlest Ascent Method: Crippen & Scheraga, 1971
Ads a term to the target function that forces the molecule up the path of
'gentlest ascent' from a local minimum to a saddle point.
Disadvantage:
One cannot be sure when the search of conformational space is sufficiently
complete. Not a priori way of knowing how many local extrema to expect on
the PES.
3. Dynamic Searching:

Real molecules find their conformational global minimum by fluctuating
about an ensemble of configurations within energetic reach (as determined
by the systems temperature).
 In principle then, if one can simulate the motions and fluctuations
experienced by a molecule, the global minimum will be sampled by this
simulation (eventually)
 MD approach: F=ma
By taking successive time steps, a time dependent trajectory can be
constructed for all the atoms representing the molecular motions.
By using available thermal energy to climb and cross conformational E
barriers, dynamics provides insight into the accessible conformational
stated of the molecule.
4.Monte Carlo:
Adiabatic Mapping:
 The simplest application of energy minimization methods to problems
of macromolecular dynamics has been the estimation of energy costs for
transitions between two or more conformations of a given molecule.
 If a suitable reaction coordinate can be defined to describe the
variation of conformation from one state to another, then the dynamics of
the system frequently follows directly from the dependence of
conformational energy on the reaction coordinate.
Example: Rotation of a single tyrosine ring through an angle of 180

A second class of problems that can be studied by comparing the relative
energies of different conformations is the large scale bending of proteins
and nucleic acids. If the effective energy profile along the reaction
coordinate is nearly quadratic, the restoring force can be calculated and
motion along this coordinate can be examined by diffusive models that use
the Langevin equation.
<
>
Example: Hinge bending motion of lysozyme
McCammon, Gelin, Karplus*Wolynes, 1976
<  Definition: The calculation of effective energy costs for transitions
by displacements along a reaction coordinate followed by relaxation.
 This name derives from an assumption implicit in the use of such
energies in discussions of dynamics, namely that he relaxation produced by
energy minimization would occur rapidly (adiabatically) on the time scale
of the transition during actual dynamics. The resulting energies are
perhaps best thought of as rough approximations to the enthalpic component
of the potential of mean force.
 The potential of mean force is a function of the reaction coordinate;
for each value of this coordinate, it specifies the free energy of the
system corresponding to a thermal average over all coordinates other than
the reaction coordinate. If the motion along the other coordinates is
faster than that along the reaction coordinate, the potential of mean force
is the true effective energy for the transition because the systematic
force along the reaction coordinate in the fluctuating molecule is the
negative gradient of this potential.

Procedure:
The molecular model is relaxed by extensive minimization to produce a
starting conformation. It is then deformed in a series of small steps with
energy minimization after each step to produce structures whose energies
are compared to those of the starting conformation. For large deformations
requiring many steps, problems arise if the energy minimization methods
cannot remove the accumulated stresses efficiently enough to produce
structures that are truly comparable to the starting conformation.
Normal Mode Analysis
 An extension of energy minimization methods that permits
calculation of atomic displacements and thermodynamic properties, sometimes
called harmonic dynamics.
 Once a local minimum energy structure for a molecular system has been
found, a massweighted second derivative matrix of the potential energy
function force constants is diagonalized giving #N eigenvectors and
eigenvalues. The first 6 eigenvectors have essentially zero eigenvalues
and are associated with global translation and rotation of the entire
system. The remaining 3N6 are the normal mode displacement vectors and the
associated eigenvalues are the frequencies for the internal dynamics of the
system.
 Given the normal modes and frequencies, one can compute dynamic and
thermodynamic properties. Normal mode calculations can be compared directly
with vibrational spectra, and the method is often used to adjust potential
energy function force constants.
 Provide info. only about the internal dynamics that can be well
represented by small harmonic displacements around a single local energy
minimum. However, the behavior of complex biomacromolecules may be better
described by anharmonic motion around one or more local minimums. There are
techniques that have been developed that include anharmonic corrections in
normal mode analysis(when anharmonicity is not too severe.)