INTRODUCTION

Top Abstract Introduction Theory Brownian dynamics simulations Results Discussion References

The dynamics of proteins and other biomolecules are essential for their biological function (Subbiah, 1996). The theoretical investigation of this dynamics is very challenging because biologically relevant motions occur on a wide variety of time scalesfrom tens of picoseconds for fluctuations around equilibrium positions to seconds for protein folding. For example, the widely used molecular dynamics (MD) simulations (Brooks et al., 1988) are currently limited to the time domain of nanoseconds at most. Realistic MD simulations of biomolecules in aqueous solution involve tens of thousands of atoms and become prohibitively expensive for longer times. Despite the continual rapid increase in computer power and technical advances in simulation algorithms (Zhang and Schlick, 1994; Figueirido et al., 1997), MD simulations reaching biologically relevant times cannot be expected in the near future. In addition, such simulations would generate an overwhelming amount of information, the analysis of which would present new problems of its own. These limitations emphasize the need for theoretical advances permitting the description of the long-time dynamics in biological systems by statistical mechanical methods.

We have been developing a theory for predicting and characterizing the long-time dynamics of peptides and proteins under conditions where one or a few configurations are insufficient for describing the range of conformational states accessed during the dynamics. There are numerous instances when such flexibility is important for biological function. For example, although most proteins in their native state are compact globular objects, many proteins have rather flexible portions that may have physiological significance. A recently discovered example is the structure of the prion protein PrP^C, which contains a 98-residue amino-terminal segment that exhibits features of a flexible random coil (Riek et al., 1997). This flexible coil may enable structural transitions to PrP^Sc aggregates. Even more striking is the absence of a stable secondary and tertiary structure in the native state of the protein p21 (involved in cell cycle control) when it is not bound to its ligand(s) (Kriwacki et al., 1996). Denatured proteins, on the other hand, are totally flexible and may roam through an enormous number of different conformations during their dynamical motion and during the process of protein folding. Recent experiments indicate that a compact protein state is achieved only late in the folding process (Sosnick et al., 1996; Smith and Dobson, 1996). Likewise, small peptides are generally not characterized by single native structures, and the accessibility of several different low-energy conformations may contribute to their ability to bind to several different receptors and thereby exhibit several different physiological functions. An example with biomedical significance is the amyloid -peptides involved in the plaque formation occurring with Alzheimer's disease. These peptides do not have a dominant conformation in solution and instead preferentially adopt random-coil, -helix, or -sheet (the latter promotes plaque formation) structures, depending on the temperature, solvent conditions, and pH (Barrow and Zagorski, 1991; Kirshenbaum and Daggett, 1995). These examples, as well as many others (see Discussion), demonstrate the need for understanding the dynamics of flexible biomolecules and also show that the currently available experimental and theoretical methods are insufficient to provide such an understanding. Therefore, the goal of our theory is to remedy this deficiency by 1) developing a general method for characterizing the motion of flexible biomolecules and 2) predicting this motion on the long time scales that are pertinent to biological functions.

As noted above, MD simulations are well developed for following the detailed trajectory of a biological system for short times and modest system sizes. However, when applied to a flexible chain molecule that may traverse many different conformational states during a single trajectory and that may exhibit very disparate sequences of conformational states between different trajectories, the MD methods are clearly limited in their descriptive powers. A statistical treatment of the dynamics therefore appears more appropriate. Statistical mechanics provides a complementary method to MD simulations for probing long-time dynamics. Rather than following a large number of individual trajectories (or a very long single trajectory) and averaging over their (its) behavior, the alternative method, from the outset, studies the dynamics as averaged over a whole equilibrium ensemble of systems. Such a statistical description is provided by various time correlation functions, which may be obtained, in principle, from averaging one long or several shorter MD trajectories but which emerge naturally from a statistical mechanical theory. Examples include the P₂(t) time correlation functions measured by NMR and fluorescence depolarization methods. However, these methods provide information only on the time scale for local bond and transition dipole reorientations, respectively. Thus, they do not directly probe the slower more global relative motions of different groups, e.g., the two domains forming the scissors motion in immunoglobulins. Hence, we recognize the need for assessing which correlation functions can characterize the slow large-scale global motion of flexible denatured proteins, the relative motions of protein domains, and the flexible motions of peptides, as well as the need for developing a theory to predict these motions. A predictive analytical theory is distinguished from MD simulations, which merely implement Newton's equations numerically. Because our goal lies in describing motions on time scales that are long compared with those accessible to MD simulations, it is clear that a predictive theory must contain fundamental conceptual advances in the statistical description of long-time dynamics.

Our theory extends to flexible chain molecules the projection operator techniques developed by Mori and Zwanzig (Mori, 1965; Zwanzig, 1961), which provide a formally exact reduced description of the dynamics of many-body systems. Zwanzig and Mori show that by identifying a suitable set of slow variables the original full Liouville equation of motion can be rewritten in the form of a generalized Langevin equation (GLE), which contains explicitly only the small set of N slow variables r(t),

(1)

with f_r(t) the random forces, ₀ a frequency matrix, and () the memory function matrix. The apparent simplicity of the GLE is deceiving; the complexity of the many-body problem is hidden in the properties of the memory functions and random forces. Nevertheless, when a clear time scale separation exists between the slow and fast degrees of freedom, the solution of the GLE can be simplified considerably by invoking a Markovian approximation for the memory functions, i.e., that the memory functions decay very rapidly in time and therefore ()  (t), where  is a relaxation matrix, and (t) is the Dirac delta function. This has allowed the successful application of the GLE in many areas of condensed-matter physics, such as the dynamics of simple fluids.

The application of the GLE to the description of the dynamics of flexible chain molecules is complicated by the fact that no clear separation of time scales exists for these systems. The lack of a time scale separation, combined with the connectivity and flexibility of the polymer chain, leads to a whole set of challenging problems. The first problem is the choice of the slow variables themselves. A second problem is the treatment of polymer-solvent interactions, a problem that can be tackled by introducing a hydrodynamic model for the solvent (Metiu et al., 1977) and, perhaps, by including an explicit shell of solvent molecules in the set of slow variables. The hydrodynamic model describes the influence of the solvent through hydrodynamic interactions and through friction coefficients for atoms or groups of atoms in the polymer chain. The final challenge arises because the GLE describing polymer dynamics is inherently multidimensional, and therefore the memory functions are N × N matrices with poorly understood properties. The theoretical treatment of these matrices is very complicated as they contain both intermolecular and intramolecular contributions.

The first approach to these problems is due to Zwanzig and Bixon (Zwanzig, 1974; Bixon and Zwanzig, 1978). They begin from the diffusion equation with hydrodynamic interactions and project onto the slow variables chosen as the linear set of bead positions in the chain. With this choice of slow variables and a Gaussian bead-spring model with preaveraging of the hydrodynamic interactions the memory function terms become zero and the GLE reduces to the optimized or generalized Rouse-Zimm (GRZ) equations of motion. The GRZ equations extend the classical Rouse-Zimm theory for polymer dynamics by including non-nearest-neighbor interactions. Perico and co-workers (Perico, 1989; Guenza and Perico, 1993) have demonstrated that the non-nearest-neighbor contributions very substantially influence both the local and global dynamics of synthetic polymers. Moreover, they show that the time correlation functions describing the dynamics can be derived exactly within the framework of the GRZ model and that the evaluation of these dynamical quantities requires as input only equilibrium information and bead friction coefficients (Perico and Guenza, 1985). However, for any realistic interaction potentials beyond the simple bead-spring model, the memory function terms no longer vanish, and due to the enormous difficulties of their calculation, these terms have been customarily neglected in the earlier treatments.

More recently, the GRZ approach has been applied to peptides containing a single tryptophane residue, glucagon, the hormone ACTH, and their fragments (Chen et al., 1987; Hu et al., 1990, 1991, 1995). The equilibrium averages required as inputs to the dynamical theory are obtained from realistic potentials, either with rotational isomeric state calculations, using the ECEPP potentials developed by Scheraga and co-workers (Li and Scheraga, 1987) or from MD simulations with the CHARMm package. The theoretical predictions are in general qualitative, and in many cases even quantitative, agreement with experimental measurements of the decay of fluorescence anisotropy of the tryptophane residue. Although this agreement is very encouraging, the neglect of the memory functions in this work, combined with the other approximations introduced in the GLE that are described above, implies that the agreement may possibly result from a fortuitous cancellation of errors. For example, it is possible to assess indirectly the significance of the memory function terms by extending in a hierarchical fashion the set of slow variables from the set of virtual bonds in the peptide to include all backbone bonds and finally the side chain bonds. Such an expansion of the slow variable set indeed shows that the memory functions exert a significant influence on the dynamics of the ACTH(6-10) fragment and that the best agreement with experiment is obtained only when all backbone and side chain bonds are included in the set of slow variables. Hence, it is obvious that additional work is required to test the significance of the individual approximations to the GLE and thereby to instill confidence in the predictive capabilities of the GRZ model and its extensions.

The complexity of these problems is best approached by separately testing each approximation in a manner that is independent of the other approximations involved. The largest and least studied of these approximations is the neglect of the memory functions in the GRZ equations. Our recent work (Perico et al., 1993; Tang et al., 1995; Kostov and Freed, 1997) has designed an approach that evaluates only the contributions of the memory functions to the dynamics predicted by the GRZ equations. The point of departure is the diffusion equation with hydrodynamic interactions neglected. The approach of Zwanzig and Bixon applied to such a diffusion equation produces a generalized Rouse (GR) model. This diffusion equation is solved alternatively in two ways: 1) numerically by Brownian dynamics (BD) simulations and 2) analytically with a theory described below. Both the theory and the simulations use the same friction coefficients and identical, realistic potential functions, and hence the theory contains no adjustable parameters. These features, as well as the absence of an approximate hydrodynamic model for the solvent in both theory and simulations, allow the unambiguous testing of the significance of the memory functions on the dynamics.

The dynamics of the polymer are described by time correlation functions (which can be obtained from the theory for arbitrary long times). The theory expresses the time correlation functions in terms of a set of equilibrium averages and the eigenvalues and eigenfunctions of the adjoint of the diffusion operator. The eigenvalues are found by expanding the eigenfunctions in a suitable basis set and by solving the resulting generalized eigenvalue problem. We have developed a procedure for the choice of the basis functions and have shown that this choice corresponds to including the contributions from the previously neglected memory functions (Tang et al., 1995; Kostov and Freed, 1997). The calculations proceed in two steps. In the first step, the necessary equilibrium averages are obtained from the BD simulations. In the second step, the generalized eigenvalue problem is solved, and the correlation functions obtained from the theory are compared with those computed from the BD simulations. The simulations are considered to provide the numerically exact correlation functions as the simulated trajectories are quite long and therefore the statistical error is small.

The selection of a suitable basis set is the most important step of both projection operator methods and our theory. If a basis is chosen as the individual atomic positions (first-order basis), the theory reduces to the generalized Rouse (GR) theory, which completely neglects the memory functions. Adding more functions in the basis set includes contributions from the memory function matrix and therefore produces a more accurate description of the internal dynamics. In the spirit of mode-coupling theory for the treatment of long-time critical dynamics (Kawazaki, 1970; Kadanoff and Swift, 1968; Kapral et al., 1976), the additional basis functions have been taken as products of the normal modes (GR modes) obtained from the first-order basis set. Thus, we refer to the theory as the mode-coupling theory, which is a more descriptive name than the initially used term matrix expansion method.

The size of the mode-coupling basis sets, however, increases dramatically with system size, as the number of possible products of the GR modes is virtually unlimited even for small peptides. For example, the automatic inclusion of all possible trilinear and pentalinear products, used initially for alkanes, yields basis sets the size of which increases as N³ and N⁵, respectively, where N is the number of atoms. This explosive growth in the number of basis functions requires the introduction of a selection procedure that includes in the basis set only those modes that contribute most to the long-time dynamics. Such a procedure was described and tested in a previous paper studying the dynamics of glycine oligomers (Kostov and Freed, 1997). A large number of high-order products of the normal modes is sorted according to their first-order relaxation times, and only a small subset consisting of the slowest decaying terms is included in the basis set. The use of such a selection procedure is motivated by physical considerations; the slowest decaying modes are anticipated to provide the most significant contributions to the long-time dynamics. A basis set selected in this manner has been demonstrated to yield quantitative agreement with the BD simulations for all global and many local bonds in triglycine and octaglycine. This basis set scales nearly linearly with system size and leads to a considerable reduction in the number of basis functions needed to provide an accurate description of the long-time internal dynamics of peptides. The new sorting procedure therefore allows the application of the theory to the long-time dynamics of larger, biologically relevant peptides.

To further test and refine the theory, in this paper we consider the internal dynamics of the neurotransmitter met-enkephalin, which is displayed in Fig. 1. Met-enkephalin is chosen to provide a very stringent test of the theory as its flexibility is combined with a considerable degree of rigidity due to its small size and the presence of two aromatic rings. Moreover, the dynamics contains the interesting motions of phenyl groups and of side chains relative to the backbone, motions that are present in general flexible peptides. In addition, met-enkephalin has been well studied both experimentally and theoretically because of its biological importance (Li and Scheraga, 1987; Vasquez et al., 1994; Montcalm et al., 1994; van der Spoel and Berendsen, 1997; Benham and Deber, 1984; Anteunis et al., 1977).

View larger version (12K): [in this window] [in a new window]   FIGURE 1   Structure of met-enkephalin (Tyr-Gly-Gly-Phe-Met). Arrows indicate the dihedral angles that may exhibit rotations.

We first briefly review the basic theory and then present the details of the Brownian dynamics simulations employed for calculating the time correlation functions and the equilibrium averages and describe the sorting procedure used to select the basis set. The results are presented and discussed in the last two sections, which provide illustrations for both global and local motions, with the latter including various phenyl group motions. The Discussion considers the biological relevance of our computations and of peptide flexibility.

	THEORY

Top Abstract Introduction Theory Brownian dynamics simulations Results Discussion References

The basic theory is described in detail elsewhere (Perico et al., 1993; Chang and Freed, 1993; Tang et al., 1995). The following summary introduces necessary notation and concepts. Let r(t) represent the 3N Cartesian coordinates of the N peptide atoms r(t) = {r₁(t), r₂(t), ... , r_N(t)}. The starting point is the diffusion equation with hydrodynamic interactions neglected:

(2)

where (r, tr₀) is the probability density that the peptide configuration is r at time t given that it is r₀ at t = 0, U(r, t) is the potential function described in the next section, _i is the friction coefficient for atom i, k_B is Boltzmann's constant, T is the absolute temperature, and  = 1/(k_BT). The equilibrium solution of Eq. 2 is the Boltzmann distribution _eq(r) = Z¹ exp (U), where Z is the partition function dr exp (U). The time evolution of any dynamical observable g(t) is given by

(3)

where is the adjoint of :

(4)

Equation 3 can be formally solved for g(t):

(5)

The application of projection operator techniques to the diffusion Eq. 2 produces the GLE of Eq. 1 (Zwanzig, 1974; Bixon and Zwanzig, 1978). The diffusion equation, however, can be solved alternatively by considering the eigenvalue problem of the operator :

(6)

where _n and _n are the complete set of eigenvalues and eigenfunctions, respectively, of . Using Eq. 6 and expanding (r, t)/_eq(r) in the complete set of eigenfunctions _n any time correlation function f(0)g(t) can be expressed as

(7)

where ··· denotes the equilibrium average:

(8)

As the eigenvalue equation cannot be solved analytically for the complicated potentials describing the interactions in biomolecules, we approximate the exact solution using a basis set expansion:

(9)

The C_in and _n are then determined from the generalized eigenvalue problem:

(10)

where S_ij = _ij is the positive definite metric matrix,  is the diagonal matrix of (approximate) eigenvalues ₁, ₂, ... , and

(11)

where the last step follows from successive integrations by parts. The normalization condition is C^TSC = I, where I is the unit matrix. The time correlation function then becomes

(12)

It is important to note that the theory contains no adjustable parameters; it uses the same potentials and friction coefficients as the Brownian dynamics simulations to be described in the next section. The only approximation lies in the use of a finite basis set in Eq. 9.

	BROWNIAN DYNAMICS SIMULATIONS

Top Abstract Introduction Theory Brownian dynamics simulations Results Discussion References

The diffusion equation is solved numerically by considering the solutions to the equivalent Langevin equation,

(13)

where F_i = U(r)/r_i is the force acting on the ith atom, and _i(t) is a Gaussian random variable with a zero mean. The fluctuation-dissipation theorem provides the variance of _i(t) as where _ij is the Kronecker delta function and (t  t') is the Dirac delta function.

The Brownian dynamics simulations are performed using the algorithm of van Gunsteren and Berendsen (van Gunsteren and Berendsen, 1982):

(14)

where t is the time step, F'_i(n) = [F_i(n)  F_i(n  1)]/t is the finite difference approximation to the time derivative of the force, and _in is a random positional displacement taken from a Gaussian distribution with zero mean and variance 2(k_BT/_i)t.

The simulations and theory use the standard GROMOS87 (van Gunsteren and Berendsen, 1987) force field, which includes the following interactions:

(15)

The local interactions are modeled by harmonic bond, angle, and improper dihedral angle terms and by specified dihedral angles. The long-range interactions include standard Lennard-Jones and electrostatic potentials. Details about the interaction potentials and the potential parameters are given in the GROMOS manual (van Gunsteren and Berendsen, 1987) and in our earlier paper (Kostov and Freed, 1997).

The dynamics of met-enkephalin at temperature T = 300 K is studied in a solvent modeled as a structureless continuum with a viscosity  of 0.84 cp. The sole effect of the solvent is in exerting a viscous damping force and a random force on each atom of the peptide. Each peptide atom is included explicitly, with the exception of the CH₃, CH₂, and CH groups, which are treated as united atom units. The atomic friction coefficients are calculated using Stokes' law as _i = 6R_i, where {R_i} are the van der Waals radii. The simulations are performed with a time step of 2 fs. The time correlation functions and equilibrium averages are calculated from 18 trajectories starting from random initial conformations. Each 15-ns data collection trajectory is preceded by energy minimization and a 3-ns equilibration. No cutoff is used for the nonbonded interactions. The total kinetic and potential energies remain constant within small fluctuations, a feature indicating the stability of the simulations.

Choice of basis set

The selection of a basis set is the most important part of the theory as it represents the only approximation that is introduced to solve the diffusion equation. We consider polynomial basis functions of the general form _i, _ij, _ijk, _ijkm, ... , where i, j, k, m, ... are integers between 1 and N  1, and {_i} can be any suitable set of functions. If all such possible functions are included in the basis set, the basis set would be complete, and the theory would produce the exact time correlation functions. However, it is impossible to use such a basis set as it contains an infinite number of functions. Therefore, it is essential to select subsets of the complete basis set that produce accurate results. The choice of basis has been extensively discussed in our previous work (Tang et al., 1995; Kostov and Freed, 1997), and therefore only a brief summary is given here for the sake of completeness.

First-order basis

Mode-coupling basis

(16)

₁t

₂t

_nt

	RESULTS

Top Abstract Introduction Theory Brownian dynamics simulations Results Discussion References

We compute the normalized interatom distance autocorrelation functions,

(17)

where l_ij = r_i  r_j, and r_i and r_j are the position vectors for arbitrary atoms i and j in the peptide. Calculations for alkanes (Tang et al., 1995) show that similar trends are obtained for the computationally more involved P₂(t) correlation functions. The correlation functions C_ij(t) are evaluated directly from the BD simulations using a time average over 18 trajectories and are displayed with solid lines in Figs. 2-6. The figures also present the same correlation functions as computed from the theory using the first-order basis set and the mode-coupling basis set generated with the sorting procedure described in the previous section. We do not present the computations for the second- and third-order basis sets as we have demonstrated previously that these basis sets 1) include a large number of functions that are inconsequential to the long-time dynamics and 2) omit contributions from other important terms. Thus, the essential contributions from these basis sets are fully incorporated by the sorting procedure.

Global motions

The global motions can be divided into several types: the first group involves only backbone atoms, the second group includes motions of the side chains relative to the backbone, whereas the third encompasses the relative motions of two side chains. Fig. 2 displays a representative illustration of the general trends for the dynamics of long-range virtual bonds in the backbone of met-enkephalin. Fig. 2 A compares the autocorrelation function for the virtual bond connecting the C₁ and C₄ atoms as determined from the simulations, the first-order basis, and the sorting procedure with 1169 and 1345 basis functions. The theory with the first-order basis set produces poor agreement with the simulations. This is expected as the first-order basis set corresponds to the GR model, which neglects the memory functions, thereby underestimating the local stiffness and predicting a too-fast dynamics. Fig. 2 A demonstrates that the basis set selected by using the sorting procedure produces very good agreement with the simulations. The quite small difference between the predictions of the theory with 1169 and 1345 basis functions indicates that convergence is practically reached and that adding more basis functions beyond 1345 would not change the results substantially. The results in Fig. 2 A, as well as those for many other correlation functions indicate that the theory reliably and consistently predicts the relaxation time scales for the long- and mid-range motions in the peptide. However, as discussed in the following subsection, the agreement between theory and simulations can further be improved by using a modified basis set that treats the rigid portions of the peptide more accurately.

View larger version (29K): [in this window] [in a new window]   FIGURE 2   Comparison of correlation functions for long-range motions in met-enkephalin as calculated directly from the BD simulations (solid line) and the first-order theory (long dashed line). (A) The C1---C4 virtual bond; comparison of the mode-coupling theory with the 1169 slowest modes (short dashed line) and the mode-coupling theory with the 1345 slowest modes (dotted line). (B) The C1---C4 virtual bond; comparison of the mode-coupling theory without centers of friction and the 1345 slowest modes (short dashed line) and the CF mode-coupling theory with the 1120 slowest modes (dotted line). (C) The end-to-end C1---C5 virtual bond. The line pattern is the same as in B.

Centers of friction (CF) basis set

A large portion of met-enkephalin consists of rigid units: the two aromatic rings as well as the peptide bonds. The local dynamics of these rigid units are the most difficult for the theory to describe (see the following subsection). In addition, the presence of rigid groups influences the global dynamics.

As the rigid units experience highly collective motions, a natural approach to treat the dynamics of these units is to combine their atoms in one group by including in the basis set only the center of friction of the group. For example, grouping the CO-NH atoms in a peptide bond gives the center of friction as

(18)

The centers of friction (R_cf^ring) of the aromatic rings are defined similarly to include all atoms in the ring along with the coplanar C united atom groups.

The bonds between neighboring atoms within the peptide groups and the aromatic rings are replaced in the first-order basis set by the bonds Cⁱ---R_cf^pep_i and R_cf^pep_i---Cⁱ⁺¹ for each peptide group and by a bond between the C atom and R_cf^ring for the aromatic amino acids. (We find that the best results are obtained when including in the basis the centers of friction of both the peptide bonds and the aromatic rings.) This reduces the total number of first-order basis functions to 20 from the original first-order basis with 39 functions. The mode-coupling basis set is then constructed as outlined above using this modified first-order basis. Spatial resolution is not lost by using centers of friction as Eq. 12 can still be used to compute the dynamics of any bond within a rigid group, even if this bond is not included explicitly in the first-order basis. Basis sets employing the centers of friction will be referred to as CF basis sets both in first-order as well as for the mode-coupling basis sets.

Fig. 2 B displays (for the same bond presented in Fig. 2 A) a comparison of the correlation functions obtained when the bonds in the rigid units are included explicitly in the first-order basis and when these units are represented by their centers of friction. Clearly, grouping together the atoms of the rigid units yields better agreement between theory and simulations. Figure 2 C presents the calculations for the end-to-end C₁---C₅ bond, which exhibit a similar trend.

The improvement in the theoretical predictions by using centers of friction is more pronounced for correlation functions corresponding to global virtual bonds between the backbone and the side chains or for side chain-side chain virtual bonds. The examples in Figs. 3 and 4 reaffirm the conclusion that the theory can predict quite accurately the relaxation times for a wide variety of long- and mid-range motions in met-enkephalin. Fig. 3 presents the calculations for two virtual bonds between the backbone and the ends of the Phe (Fig. 3 A) and Met (Fig. 3 B) side chains. Fig. 3 A emphasizes again that using the CF basis set yields a substantial improvement whereas Fig. 4 A displays the correlation functions for a virtual bond connecting the two aromatic rings, where the contributions from the CF basis are most pronounced. On the other hand, the correlation function in Fig. 4 B for the relative motions of the Phe and Met side chains incurs negligible improvement from using the CF basis. Although the CF basis set leads to improved theoretical predictions for many bonds, in some cases, as for example the bonds in the Met side chain, the CF basis yields identical correlation functions to those from the original mode-coupling basis set without the centers of friction grouping. As this lack of improvement usually appears for bonds that are not part of the rigid units, these results can generally be explained by the absence of additional coupling between the motions of the rigid units and these particular bonds.

View larger version (30K): [in this window] [in a new window]   FIGURE 3   Comparison of the correlation functions in met-enkephalin as calculated from the simulations (solid line), the first-order theory (long dashed line), the CF mode-coupling theory with 1120 modes (short dashed line), and the mode-coupling theory without centers of friction and 1345 modes (dotted line) for (A) the C1---C4 backbone-Phe bond and (B) the C1---C5 backbone-Met bond.

View larger version (31K): [in this window] [in a new window]   FIGURE 4   Comparison of the correlation functions in met-enkephalin for (A) the C1---C4 Tyr-Phe virtual bond and (B) the C4---C5 Phe-Met bond. The line pattern is the same as in Fig. 3.

Two advantages accrue from using a centers of friction basis set for the rigid groups. 1) The number of basis functions in the first-order basis is significantly reduced, thereby allowing the treatment of larger polypeptides. Although the twofold reduction from 39 to 20 functions is inconsequential for met-enkephalin as the number of the mode-coupling basis functions is far greater than the size of the first-order basis, this size reduction may become important for proteins where the first-order basis is already large. 2) The theoretical predictions are improved for the correlation functions of a large number of bonds, in particular for bonds involving the rigid units. Therefore, hereafter we present only calculations with the CF basis set using 1120 mode-coupling terms.

Local Motions

Previous work for glycine oligomers (Kostov and Freed, 1997, 1998) demonstrates that the theory provides a good representation for the long-time decay of many local correlation functions. The present work for met-enkephalin confirms these trends, as depicted in Fig. 5 A for the N₃---C₃ bond. However, the accurate prediction of the local dynamics for some of the local bonds represents a challenge for the theory. This is due, in part, to the fact that the theory is developed to treat the dynamics of flexible polymers, whereas the local rigidity of the peptide chain is manifested to a greater extent for the local than for the global motions. In addition, our method is directed at the long-time dynamics, whereas many of the local bonds (such as N---H bonds) exhibit fast relaxation. Thus, computations for a significant number of local bonds appear to exhibit a considerable discrepancy between the theory and BD simulations. The differences in the relaxation times, as obtained from integrals of the correlations functions ( = ₀ dtC(t)), vary from a factor of 2---3 for backbone bonds that are part of the peptide group to a factor of 5 or more for the side chain bonds involving the aromatic rings. Such a discrepancy, although to a lesser extent, is contained in previous computations for the oligoglycines. Its origins are better evidenced by semi-log plots of the correlation functions as presented in Fig. 5 B for the C₁---N₂ bond. Fig. 5 B clearly demonstrates that the theory, indeed, correctly predicts the long-time relaxation rates of the correlation functions (i.e., the slopes of the theoretical curves agree with the simulations) and that the discrepancy between theory and BD simulations for the correlation functions is solely due to an initial very fast decay in the theory at short times of up to a few hundred picoseconds. This initial fast decay in the theory is simply explained by noting that the selection of the slowest modes by the sorting procedure emphasizes the long-time dynamics at the expense of omitted modes responsible for accurately depicting the shorter-time dynamics.

View larger version (27K): [in this window] [in a new window]   FIGURE 5   Comparison of the correlation functions in met-enkephalin as calculated from the simulations (solid line), the first-order CF theory (long dashed line), and the CF mode-coupling theory (1120 modes) (short dashed line) for (A) the N3---C3 (3) bond. (B) Semi-log plot of the C1---N2 bond.

More significantly, the largest discrepancies between theory and simulations are seen for bonds that are part of rigid units, such as the peptide bond or the phenyl side groups. The strongly correlated motions exhibited by such rigid units are difficult to capture with the polynomial basis set. Yet, even in these cases, the long-time portions of the theoretical correlation functions correctly predict the relative relaxation times for the different motions, e.g., the motions of the aromatic rings in different directions. Correlation functions for a set of bonds in the phenyl rings are displayed in Fig. 6, which reaffirms the assertion that agreement between theory and simulations is obtained for the long-time relaxation rates.

View larger version (28K): [in this window] [in a new window]   FIGURE 6   Semi-log plot comparison of the correlation functions in met-enkephalin for (A) the C1---C1 bond in the Tyr ring and (B) a C---C bond in the Phe ring. The line pattern is the same as in Fig. 5.

The short-time deficiency of the theory is natural and is not of great concern as the short time scales are readily accessible by molecular dynamics simulations. Nevertheless, we have tried to improve the results at short times by selecting alternative basis sets. One possible approach is to note that the motion of a particular local bond of interest is most likely to be correlated to the motions of its neighboring bonds. We can therefore form a basis set focused on describing the dynamics of a specific bond by including in the basis set products of the GR modes (describing the global motions) together with products of the bond vectors for the chosen bond and its neighbors (describing the local correlations), where the Greek indices refer to the selected bond and its neighbors. The addition of the neighboring bond correlations indeed produces some improvement in the correlation function for the selected bond at short times. This improvement, however, is rather small and does not alter substantially the agreement with the BD simulations.

A similar approach, termed the maximal correlation mode-coupling approximation (MCA), has been recently proposed by Perico and Pratolongo (1997). They study as model systems the dynamics of the freely jointed chain, broken rods, and rods and show that the dominant corrections to the optimized Rouse-Zimm dynamics (a first-order theory) are contained in a small number of product functions involving only the bonds having maximal correlation with the bond of interest. Their MCA basis set adds to the first-order basis all diagonal trilinear products of the form and thus is twice as large as the first-order basis set. A variation of this basis set for treating the dynamics of a given local bond l involves including all trilinear nondiagonal terms containing that bond. Perico and Pratolongo show that for the rods this extended basis set improves the MCA predictions by only a few percent. Although the MCA basis contains the dominant corrections to the first-order basis set for the freely jointed chain and rods, our calculations demonstrate that for met-enkephalin this MCA basis improves the first-order dynamics to a lesser extent than achieved by adding to the first-order basis the same number of normal mode products selected with the sorting procedure. In addition, when applied to more complicated peptides and polymers, the MCA approach has the drawback that a different basis set must be constructed for every bond for which the dynamics is studied. Therefore, even though the dimensionality of the basis set is reduced, the generalized eigenvalue problem must be solved many times if the dynamics are of interest for a large number of bonds. In contrast, our mode-coupling basis set provides a universal description of the dynamics of all global and local bonds in terms of a single set of basis functions (and eigenvectors).

	DISCUSSION

Top Abstract Introduction Theory Brownian dynamics simulations Results Discussion References

The global and backbone dynamics of met-enkephalin exhibit many similarities when compared with the previously studied dynamics of oligoglycines (Kostov and Freed, 1997). The relaxation time scales in met-enkephalin are longer by ~50% as compared with octaglycine. Although quantitative agreement of the theory with simulations is obtained for the oligoglycines, the calculations for met-enkephalin display a small discrepancy with the simulations. Semi-log plots indicate that this discrepancy is almost entirely due to an incorrect initial fast decay in the theory. This discrepancy is specific to bonds in rigid groups and might arise from the particular structure of met-enkephalin, which has a large number of rigid units relative to the size of the peptide. Future studies are necessary to determine whether this discrepancy is present for peptides with other lengths and sequences. Met-enkephalin has many interesting motions of phenyl groups and of the side chains relative to the backbone that are absent in the simpler oligoglycines.

Additional insight into the internal motions of met-enkephalin may be obtained by comparing the slopes of the correlation functions at long times. The slopes for various correlation functions differ, frequently by as much as 50%. This indicates that the overall rotation of the molecule, expected to provide the dominant long-time contribution, does not solely determine the very long-time dynamics of both the global and the local bonds. Rather, the overall rotational motion is coupled to other motions on shorter length scales in a manner that is specific for each bond. The local motions in met-enkephalin are, therefore, not only determined by their local environment, but also strongly coupled to more global motions on larger scales. In contrast, traditional interpretations of NMR relaxation data, as for example the model of Lipari and Szabo (1982), frequently rest on the assumption that the local motions in polypeptides are decoupled from the global rotational motion of the molecule. Although this assumption may be appropriate for proteins, the present computations show that it is not valid for small peptides, and additional studies of this assumption are warranted for very flexible portions of proteins. Our theoretical prediction of the coupling between global and local motions accords with NMR experiments for triglycine (Daragan and Mayo, 1994).

The theoretical correlation functions presented in Figs. 2-6 are typically obtained with more than 1000 basis functions. Therefore, the correlation functions are in general represented by a sum of more than 1000 exponentials, although many amplitudes are vanishingly small and only a few terms contribute at the longer times. This number is substantially greater than the two or three exponentials typically used to fit experimental data. However, the correlation functions computed with the theory do not represent a fit. As stressed several times, the theory contains absolutely no adjustable parameters. The same set of eigenvalues present in the exponents for the corr