Molecular Mechanics Minimization:

• Minimization
• Exploring Multidimensional PES
  1. MD:
  2. MC
  3. 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?
  
  1. Comparison of structures/properties
  2. 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.
  3. Template forcing
  4. Systematic mapping of E space
  5. Binding energies
  6. Docking
  7. Harmonic analysis
  8. 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
  • Conjugate-direction
  • Newton methods
    1. Quasi-Newton-Raphson
    2. Full Newton-Raphson
  • Quenched dynamics
  
  Minimum information to start minimization:
  1. Starting set of atomic coordinates
  2. Parameters for various terms of the potential energy function
  3. Description of molecular topology (covalent structure) (potential energy function)
  4. Additional features
     1. Use of periodic boundary conditions so that the simulated system is effectively part of an infinite system.
     2. Forces of constraint to hold the system in a particular configuration or move it toward a different target conformation.
     3. 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:
  1. Assist in refinement:
     • Structures of peptides (Scheraga, 1968)
     • Crystal structure of proteins (Levitt and Lifson, 1969; Levitt, 1974)
     • Nucleic Acids (Jack et al., 1976)
  2. Used to optimize bond lengths and angles in crystallographic refinement procedures. (Jack and Levitt, 1978; Konnert and Hendrickson, 1980)
  3. 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).
     
     
  4. 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 non-bond 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: Protein-substrate-solvent system with 5000 atoms (typical of a small protein (100-150 residues) surrounded by a 1-2 ;layers of water. Figure shows the number of nonbonded pair-wise 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 [Ala-Pro-D-Phe]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:
      1. Rigid body minimization
      2. Torsional minimization
      3. 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 4-14 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.
      
      1. 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
      2. 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
      3. 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:
      1. More control over the minimization
      2. Exploring defined regions of the conformational space
      3. Testing conformational hypothesis (template forcing;Brooks et al.,1983)
      Types:
      1. Torsional Restraints
      2. Distance Restraints
      3. Template Forcing
      4. Tethering
      
      1.Torsonial Restraints:
      • Places a harmonic torque about a bond to force it to a new value
        
        Etorque=ktorque(0-0target)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(Rij-Rtarget)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 NOE-type 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:
      1. Locate 'local', not 'global' minimizers.
      2. Designed to 'ignore' configurations if E increases
      3. 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:
        1. Screen from a large number of systematically generated points in conformational space, for the 'best' starting points from which to begin a minimization.
        2. Modify target function and/or minimization algorithms to allow systems to overcome local barriers enroute to finding lower energy minima.
        
        
        Methods:
        1. Systematic searching
        2. "Uphill minimizers
        3. Minimization through reduction of dimensionality
        4. Dynamic searching
        5. 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:
        1. Crippen & Scheraga, 1971; 1973 Employ pruning algorithms that ignore large volumes of unlikely conformational space.
        2. 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:
      1. Overcomes the natural tendency of minimizers to get stuck in local wells.
      2. Methods:
         1. 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.
         2. 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 mass-weighted 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 3N-6 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.)