Dynamics - General Methodology

Contents

• Main Uses of Dynamics Simulations
• Exploring Conformational Space using Dynamics
  • Conformational Search of ANF using Dynamics and Minimization
  • Characterizing Vasopressin Dynamics

Main Uses of Dynamics Simulations

• Conformational searches.
• Generating statistical ensembles.
• Studying the motions of molecules.

Studies of molecular motions can be used to derive properties such as diffusion coefficients.

Exploring Conformational Space using Dynamics

Real molecules find their conformational global minimum by fluctuating about an ensemble of configurations within energetic reach (as determined by the system temperature). In principle then, if the motions and fluctuations experienced by a molecule are simulated, the global minimum is sampled by this simulation (eventually). This is the approach taken by molecular dynamics. In molecular dynamics, the fact that the negative of the gradient of the potential energy is equivalent to the force is used. Given the force and mass for each atom, Newton's equation of motion (F = ma) can be numerically integrated to predict where atoms will move within a short time interval. By using successive time steps, a time-dependent trajectory can be constructed for all atoms which thus represents the molecular motions.

By using the available thermal energy to climb and cross conformational energy barriers, dynamics provides insight into the accessible conformational states of the molecule. Molecular dynamics has been used in a variety of small peptide systems to find lower-energy conformational states across energy barriers that would be inaccessible to classical minimization. In particular, it has been applied to the study of vasopressin (Hagler 1985) and gonadotropin-releasing hormone (GnRH) (Struthers et al. 1984, Hagler 1985).

For GnRH, novel constrained analogs with 100-fold higher potency than the parent cyclic compounds were designed, based on the conformations of GnRH and GnRH analogs identified with a combined strategy incorporating energy minimization, template forcing, and molecular dynamics.

One of the disadvantages of dynamics, however, is that many hundreds and even thousands of conformational states are sampled during a trajectory, thus compounding the analysis process. Furthermore, the transitions of interest may take place on a time scale of milliseconds to minutes or longer at 300° and, at present, dynamic simulations are limited to tens of nanoseconds. Minimizing periodically during a dynamics trajectory helps the analysis problem by identifying those structures that represent the local minimum about which the dynamics fluctuations take place. This is illustrated in Figure 5-7, where the dark magenta line superimposed on a simple elliptical energy surface represents a dynamics trajectory of a particle beginning at rest from point a. The initial gradient pulls the particle downhill but, because energy is conserved, its velocity increases as it gets closer to the minimum. Its momentum takes the particle back and forth, visiting many points on the energy surface during the trajectory. However, if we stop the particle at any point during the dynamics and minimize, only one point is found. In higher-dimensional systems, such as real molecules, minimization from various dynamics points may result in the identification of several different local minima. However, the number of minima is far fewer than the number of points sampled during dynamics, and these minima are useful as points of reference in analyzing molecular structures and energetics.

Of course, a real molecular system has far too many degrees of freedom to allow plotting the trajectory on a realistic energy surface. Any attempt to project a trajectory into one or two dimensions results in a far more complex pattern that obscures which points, if any, are close to a minimum. In real calculations, then, periodic (systematic) minimizations must be used to discover which minima are being sampled by dynamics. A detailed discussion of how dynamic searching has been used to characterize the conformational states of the peptide hormone, atrial natriuretic factor (ANF) is presented as an illustration.

Conformational Search of ANF using Dynamics and Minimization

The strategy included these steps:

1. System setup
2. Classification of conformational states
3. Refinement of strained structures via annealing
4. Applying simulation results to specific design strategies

Mackay et al. (1989) performed a set of molecular dynamics calculations on ANF at elevated temperature to search for structural features that may be significant for ANF binding to its receptor. The goals of these calculations were: (1) to search the conformational space of the molecule; (2) to find unique regions that may be cross-linked to test conformational hypotheses; and (3) to compare different approaches using dynamics to search conformational space.

1. System setup and strategic approach

   The principal active form of ANF in the blood is a 28-amino acid fragment with a disulfide bond between residues 7 and 23:

                    5                   10
   SER-LEU-ARG-ARG-SER-SER-CYS-PHE-GLY-GLY-ARG-MET-
                            /                    ASP-
                           S                      ARG-
                          S                        ILE-  15
                         /                        GLY-
   TYR-ARG-PHE-SER-ASN-CYS-GLY-LEU-GLY-SER-GLN-ALA-
                25                  20
   

   An initial geometry was generated by first building the sequence into an extended conformation using standard amino acid bond lengths, angles, and side chain dihedral angles. The disulfide bridge was formed by imposing a -turn at GLY 16 and minimizing the resulting structure as a starting conformation for use in dynamics.

   Solvent was not included in this system, because the interest was not in the conformation of ANF in solution, but rather in ANF bound to its receptor. There was no explicit receptor model with which to interact; the objective was to look for any conformations accessible to the ANF regardless of the environment. High-energy structures were not ignored because the receptor may stabilize them, and the structures found were compared with those of other analogs. Common structural features among ANF and ANF analogs may suggest motifs that can be used to construct a putative active conformation.

   Since solvent was not included, all the amino acids were defined to be in their neutral states, to prevent unscreened interactions between charged side chains from dominating the structural search. Strong salt bridges often form that constrain the molecule from finding alternative conformations during dynamics. This structure was then equilibrated by running dynamics at 300 K and 600 K for 20 ps each. Data were collected from a subsequent 100 ps run at 900 K. A total of 100,000 configurations (one every 10^-15 seconds) was sampled during this simulation. To reduce the volume of data to a more manageable level, instantaneous dynamic structures were minimized at one-picosecond intervals, reducing the number of structures to be analyzed to 100.

2. Classification of conformational states

   Similarity of the structures minimized along the trajectory was quantified by the rms deviation between backbone atoms of residues 7-23 (the ring residues). The similarity between structures was evaluated by comparing the rms deviations between each possible pair of structures. Then the rms deviations were mapped in a gray scale and plotted in a 100 x 100 matrix. Elements near the diagonal represent structures that are close in time, which generally had low rms values. Off-diagonal elements represent structures separated by long times, which generally had high rms values. A number of boxlike areas appeared along the diagonal, indicating that several distinct families of structures were encountered during the progress of the dynamics calculations. The 100 structures found during dynamics were classified into 11 conformational groups that had mutual rms deviations less than a threshold value (3 Å) and were located sequentially in time.

   High temperatures greatly increase the efficiency of producing conformational transitions. Consider a simple Arrhenius model for reaction rates:

   Eq. 5-30:

   

   where k is the reaction rate, A is a pre-exponential factor describing the frequency of thermal collisions, E is the activation energy, R is the gas constant, and T is the temperature. Assuming a value for A of 10¹² sec^-1, which is a typical rate of collisions for liquid rare gas atoms, the E that yields a reaction rate of once per 100 ps at 300 K is about 5 kcal mol^-1. At 900 K, this barrier is crossed 250 times in 100 ps.

   The conformational search of ANF at 900 K produced at least 11 distinct families in 100 ps. A similar control run at 300 K produced no significant new conformational states over the same period. Dynamics at 900 K was not without risk, however. For example, the peptide bond was observed to have undergone trans-cis conversions in 42 of 100 minimized structures. This is clearly a result of the high temperature--no trans-cis transitions were observed in room temperature simulations. To avoid formation of cis bonds, a torsional restraint of 5 kcal mol^-1 rad^-2 was added to the peptide bond, which effectively eliminated trans-cis conversions. This is another good example of how restraints can be incorporated into modeling strategies. Without a restraint to prevent trans-cis interconversion, high-temperature dynamics would not be as useful for conformational searching of peptides.

3. Refinement of strained structures via annealing

   A second side effect of high-temperature dynamics is that it leads to higher-energy minima--starting from a higher point on the energy surface results in a higher-energy local minimum. By subsequently subjecting the minimized structures to lower-temperature dynamics at 300 K for 10-20 ps, up to 60 kcal mol^-1 of strain could be removed. This annealing process also changed the relative order of the 11 conformational families and was therefore crucial for an energy-based selection of conformations.

   Overall, minimization and structural classification strategies were used to refine 100,000 dynamic structures to 11 representative structures. This more manageable set of conformations was then analyzed in detail for important structural features to be incorporated into a strategy for designing analogs.

4. Applying simulation results to specific design strategies

   The biological significance of a structure can best be assessed through structure-activity relationships. Typically, the presence and location of specific functional groups are correlated with either decreased or enhanced activity. This correlation can be used to suggest requirements for activity. However, constraining conformational features should also be effective in distinguishing the biological significance of conformational families. Conformational features common to active analogs but inaccessible to inactive ones may be required for activity. Introducing a chemical cross-link between two residues can constrain the resulting structure to remain in this conformation. Identifying a unique residue-residue interaction for each of the 11 identified ANF families and introducing a chemical cross-link to restrain this residue-residue distance could provide new analogs for use in testing such conformational features. A matrix of C-C distances was calculated for each of the eleven representative annealed structures. These unique distances might be used as the basis for such a strategy (Rodgers, to be submitted).

Characterizing Vasopressin Dynamics

The strategy includes these steps:

1. System setup and strategy
2. Structural and energetic analysis
3. Mapping an adiabatic torsion potential using torsion restraints

Thus, its structure includes a constrained ring of six residues and a flexible three-residue tail. Hagler et al. (1985) reported a dynamics study of vasopressin with the objective of characterizing the dynamic transitions and conformational equilibria of this flexible peptide.

As in the ANF example, minimization techniques play a major role in understanding and analyzing dynamic results.

1. System setup and strategy

   The dynamics of vasopressin were carried out in vacuo and with neutral side chains for the same reasons as for ANF. The initial structure was modeled after a postulated structure of the homologous peptide oxytocin. Dynamics was run at room temperature for approximately 8 ps, and the resulting trajectory analyzed for transitions between conformational states. Similar to the ANF study, minimization was performed periodically during the trajectory to characterize the conformational states being sampled.

2. Structural and energetic analysis

   The flexibility of vasopressin was immediately confirmed by the dynamics. Surprisingly, even residues within the six-residue ring underwent several backbone transitions during the short simulation. Phe³, for example, started in a C₇^eq state and proceeded to a right-handed -helix before eventually finding a C₇^ax state, all in just over 6 ps. As was done for ANF, structures found during the dynamics run were minimized to confirm that the dynamics trajectory was passing through legitimate minima and not just fluctuating about a single minimum. The instantaneous conformational states, along with the resulting minimized states, are shown in Table 5-2 for Phe³. Note that dynamic structures as characterized by the (, ) values for Phe³ can vary widely (by as much as 100 degrees) but still minimize to the same conformational state. Nevertheless, there are three distinct states to which the dynamic structures minimize. This rigorously demonstrates that transitions are occurring.

   One of the most significant results of this study was that the energy levels for conformational states of Phe³ in vasopressin are more closely spaced than for isolated phenylalanine. For example, the difference between C₇^eq and the -helical state in an isolated blocked Phe residue is calculated to be 4.7 kcal mol^-1, while in vasopressin it is only 2.2 kcal mol^-1. This effective leveling of the energy levels imposed by constraints in the ring contributes to an enhanced flexibility in this part of the vasopressin molecule.

3. Mapping an adiabatic torsion potential using torsion restraints

   As a final example of how minimization can be used in conjunction with restraints to perform energy analysis, consider once again the rigid geometry versus flexible geometry torsion map. As pointed out above, energy barriers to rotation calculated by rotating about a bond with rigid bonds and angles can be misleading. What would appear to be very high barriers due to van der Waals overlap are much lower when small bond stretches and angle bends can absorb the strain.

   To illustrate this point, consider Figure 5-9, which superimposes the Phe³ (, ) trajectory onto an adiabatic (that is, flexible-geometry) contour surface. Note that the barrier, particularly from the -helix to the C₇^ax conformations, is only 3-4 kcal mol^-1. The flexible-geometry torsion map is particularly useful for understanding a dynamics trajectory, because it can nicely explain the path the molecule follows on the energy surface.

   The Constraints and Restraints section includes torsion restraints, which are used to prepare adiabatic torsional potential energy surfaces such as the , map used in Figure 5-9.

Main access page

Theory/Methodology access

Dynamics access Dynamics - Constraints References