Normal Modes and Essential Dynamics

1.1 Standard Normal Mode Analysis

Normal mode analysis (NMA) is one of the major simulation techniques used to probe the large-scale, shape-changing motions in biological molecules [1–3]. Although it has connection to the experimental techniques of infrared and Raman spectroscopy, its recent application has been to predict functional motions in proteins or other biological molecules. Functional motions are those that relate to function and are often the consequence of binding other molecules. In NMA studies, it is always assumed that the normal modes with the largest fluctuation (lowest frequency modes) are the ones that are functionally relevant, because, like function, they exist by evolutionary design rather than by chance. The ultimate justification for this assumption must come from comparisons with experimental data and indeed studies that compare predictions of an NMA with transitions derived from multiple x-ray conformers do suggest that the low-frequency normal modes are often functionally relevant.

NMA is a harmonic analysis. In its purest form, it uses exactly the same force fields as used in molecular dynamics simulations. In that sense, it is accurate. However, the underlying assumption that the conformational energy surface at an energy minimum can be approximated by a parabola over the range of thermal fluctuations is known not to be correct at physiological temperatures. There exists abundant evidence, both experimental [4] and computational [5], that the harmonic approximation breaks down spectacularly for proteins at physiological temperatures, where, far from performing harmonic motion in a single energy minimum, the state point visits multiple minima crossing energy barriers of various heights. Thus, when performing NMA, one has to be aware of this assumption and its limitations at functioning temperatures.

A standard NMA requires a set of coordinates, a force field describing the interactions between constituent atoms, and software to perform the required calculations. The performance of an NMA in Cartesian coordinate space requires three main calculation steps: 1) minimization of the conformational potential energy as a function of the atomic Cartesian coordinates; 2) the calculation of the so-called “Hessian” matrix, which is the matrix of second derivatives of the potential energy with respect to the mass-weighted atomic coordinates; and 3) the diagonalization of the Hessian matrix. This final step yields eigenvalues and eigenvectors (the “normal modes”). Each of these three steps can be computationally demanding, depending on the size of the molecule. Usually, the first and final steps are the bottlenecks. Normally, energy minimization is demanding of CPU time and diagonalization is demanding of CPU time and memory because it involves the diagonalization of a 3N × 3N matrix, where N is the number of atoms in the molecule. We have called this NMA “standard” NMA to distinguish it from the elastic network model NMA.

1.2 Elastic Network Models

Because of the computational difficulties of standard NMA, the current popularity of the elastic network models is not surprising. This is still an NMA, but the protein model is dramatically simplified. Tirion first introduced it into protein research [6]. As the name suggests, the atoms are connected by a network of elastic connections. The method has two main advantages over the standard NMA. The first is that there is no need for energy minimization because the distances of all of the elastic connections are taken to be at their minimum energy length. Second, the diagonalization task is greatly reduced compared with the standard NMA method because the number of atoms is reduced from the total number of atoms to the number of residues, if one uses only C^α atoms, as is common practice. This leads to a tenfold reduction in the number of atoms. Unlike standard NMA, elastic network models have two parameters to be set. One is the force or spring constant, normally denoted as γ or C, and the other is a cut-off distance, denoted R_c.

A pertinent question is whether the method is any less accurate than the standard NMA. Tirion showed that there is a respectable degree of correspondence between the two methods [6]. Given the drastic assumptions that are inherent in the standard NMA, the small difference between the results from these two methods is probably unimportant relative to differences between standard NMA and reality. Comparisons between movements in the low-frequency modes derived from elastic network models of 20 proteins with movements derived from pairs of x-ray structures [7] suggest the same level of moderate correspondence seen in similar studies using standard NMA. This, together with the relatively low computational cost of elastic network models, explains their current popularity in comparison with standard NMA.

1.3 Essential Dynamics and Principal Components Analysis

Because of the complexity of biomolecular systems, molecular dynamics simulations can be notoriously hard to analyze, rendering it difficult to grasp the motions of interest, or to uncover functional mechanisms. A principal components analysis (PCA) [8–10] often alleviates this problem. Similar to NMA, PCA rests on the assumption that the major collective modes of fluctuation dominate the functional dynamics. Interestingly, it has been found that the vast majority of protein dynamics can be described by a surprisingly low number of collective degrees of freedom [9]. For the analysis of protein molecular dynamics simulations, this approach has the advantage that the dynamics along the individual modes can be inspected and visualized separately, thereby allowing one to filter the main modes of collective motion from more local fluctuations. Because these principal modes of motion could, in many cases, be linked to protein function, the dynamics in the low-dimensional subspace spanned by these modes was termed “essential dynamics” [9], to reflect the notion that these are the modes essential for function. The subspace spanned by the major modes of collective fluctuations is accordingly often referred to as “essential subspace.” The fact that only a small subset of the total number of degrees of freedom dominates the molecular dynamics of biomolecules not only aids the analysis and interpretation of molecular dynamics trajectories, but also opens the way to enhanced sampling algorithms that search the essential subspace in either a systematic or exploratory fashion [11–14].

In contrast to NMA, PCA of a molecular dynamics simulation trajectory does not rest on the assumption of a harmonic potential. In fact, PCA can be used to study the degree of anharmonicity in the molecular dynamics of a simulated system. For proteins, it was shown that, at physiological temperatures, especially the major modes of collective fluctuation are dominated by anharmonic fluctuations [9, 15]. Overall, protein dynamics at physiological temperatures has been described as diffusion among multiple minima [16–18]; on short timescales, the dynamics are dominated by fluctuations within a local minimum (that can be approximated well by a system's local normal modes), whereas, on longer timescales, the large fluctuations are dominated by a largely anharmonic diffusion between multiple wells.

In NMA the modes of greatest fluctuation are those with the lowest frequencies. As in PCA, no assumptions are implied regarding the harmonicity of the motion, modes are usually sorted according to variance rather than frequency. Nevertheless, the largest-amplitude modes of a PCA usually also represent the slowest dynamical transitions.

2.1 Standard NMA

NMA is usually performed in a vacuum, where the potential energy of a biomole-cule is a complex function of its 3N coordinates, N being the number of atoms. This function is normally written in terms of its bonded and nonbonded energy terms. It is usual to use Cartesian coordinates [3], although dihedral angles have been used [1,19,20]. The basic idea is that, at a minimum, the potential energy function V can be expanded in a Taylor series in terms of the mass-weighted coordinates q_i = √m_i Δx_i, where Δx_i is the displacement of the ith coordinate from the energy minimum and m_i is that mass of the corresponding atom. If the expansion is terminated at the quadratic level, then because the linear term is zero at an energy minimum:

Thus, the energy surface is approximated by a parabola characterized by the second derivatives evaluated at the energy at the minimum. The basic, but false, assumption of NMA of biomolecules at physiological temperatures is that fluctuations still occur within this parabolic energy surface. It is known, however, that, at these temperatures, the state point moves on a complex energy surface with multiple minima, crossing energy barriers of various heights [4]. The second derivatives in Eq.1 can be written in a matrix, which is often called the “Hessian,” F. Determination of its eigenvalues and eigenvectors (equivalent to diagonalization) implies: where w _j is the jth eigenvector and ω² _j is the jth eigenvalue. There are 3N such eigenvector equations. Each eigenvector specifies a normal mode coordinate through:

The sum is over the elements of w _j. Note that ǀw _jǀ = 1. It can be shown that these normal mode coordinates oscillate harmonically and independently of each other each with the angular frequency, ω_j:

Here, A_j is the amplitude and ε_j is the phase. These normal mode coordinates are collective variables because they are linear combinations of the atom-based Cartesian coordinates, as shown in Eq.3. If a single normal mode j is activated, then: which means that, in the jth mode, the relative displacements of the Cartesian coordinates are specified by the elements of w _j. Each normal mode then specifies a pattern of atomic displacement. For example, in a multidomain protein, this pattern of displacement could indicate the relative movement of two domains. Figure 1b shows an example. A more thorough introduction to the theory and its application to biomolecules can be found elsewhere [21].

It can be shown that the lower the frequency, the larger the fluctuation of the corresponding normal mode coordinate [22]. It is common to compare the lowest frequency modes with functional modes derived from, e.g., a pair of x-ray structures, one bound to a functional ligand and the other unbound. The overlap with the jth mode can be defined as [23]:

2.2 Elastic Network Models

There is, in fact, no essential difference between the elastic network NMA and the standard NMA other than the force field. In the case of the elastic network, the Hessian would be derived from the following potential energy function [6]: where r_ij is the distance between atoms i and j and r⁰ _ij is the distance between the atoms in the reference structure, e.g., the crystallographic structure. This summation is only performed over atoms less than a cut-off distance R_c, and γ is the spring, or force constant for the elastic bond between the atoms and is the same for all atoms pairs (see Fig. 1a). The energy function in Eq.7 seems to be the most popular, although other types of functions can be used. The network corresponding to the energy function of Eq.7 is sometimes referred to as the anharmonic network model [24]. A Gaussian network model has a different energy function, which results in modes without any directional information [24,25] and will not be considered here. Once the function of Eq.7 has been calculated, the procedure is exactly the same as for the standard NMA, namely, the Hessian is calculated and its eigenvalues and eigenvectors are determined. Whereas the standard NMA must be performed on all atoms as required by the force field, the elastic network model can be carried out on a subset of atoms. Often, for a protein, this would be the C^α atoms. Compared with the standard NMA, this would result in a Hessian approximately tenfold lower in order, thus, yielding considerable computational savings in the calculation of the eigenvalues and eigenvectors because these routines are normally of the order of N³ operations, where N is the order of the Hessian matrix.

2.3 Essential Dynamics and PCA

After superposition to a common reference structure, a variance—covariance matrix of positional fluctuations is constructed: where <> denotes an ensemble average. The coordinates x are denoted as a function of time for clarity, but may be provided in any order and can be, for example, a molecular dynamics trajectory or a set of experimental structures. C is a symmetric matrix that can be diagonalized by an orthogonal coordinate transformation T: with Λ the diagonal (eigenvalue) matrix and T containing, as columns, the eigenvectors of C. The eigenvalues λ, correspond to the mean square eigenvector coordinate fluctuation, and, therefore, contain the contribution of each principal component to the total fluctuation. The eigenvectors are usually sorted such that their eigenvalues are in decreasing order. For a system of N atoms, C is a 3N = 3N matrix. If at least 3N configurations are used to construct C, then 3N – 6 eigenvectors with nonzero eigenvalues will be obtained. Six eigenvalues should be exactly zero, of which the corresponding eigenvectors describe the overall rotation and translation (that is eliminated by the superposition). If only M configurations are available (with M < 3N), then at most M – 1 nonzero eigenvalues with corresponding eigenvectors will result. If μ_i is the ith eigenvector of C (the ith column of T), then the original configurations can be projected onto each of the principal components to yield the principal coordinates p_i(t) as follows:

Note that the variance < pi² > equals the eigenvalue λ_i. These projections can be easily transformed back to Cartesian coordinates for visualization purposes as follows:

Two sets of eigenvectors μ and ν can be compared with each other by taking inner products:

Subspace overlaps are often calculated as summed squared inner products: expressing how much of the n-dimensional subspace of set μ is contained within the m-dimensional subspace of set ν. Note that m should be larger than n to achieve full overlap (O = 1).

3.1 Standard NMA

For standard NMA, one needs a set of coordinates, a force field, and software to perform the calculations. Often NMA is performed using molecular mechanics software packages that are also able to perform molecular dynamics simulations, etc. For a protein, the structural information is normally held in a PDB file. The software will normally be able to interpret the file to determine the correct energy function using the selected force field. Any missing atoms should be added. Missing hydrogen atoms also need to be added but most software packages have routines to do this. It is usual, but not a requirement of the methodology, to remove water and ligands. Once the system is prepared, the first major calculation is energy minimization.

3.2 Elastic Network Models

From the last step onward, there is no essential difference to the standard NMA. However, if calculations are performed on the C^α atoms only, then, naturally, one can only compare with the movements of the C^α atoms in the experimentally determined functional mode, i.e., movements of side chains cannot be compared. Tama and Sanejouand have made a comparison between the results from an elastic network model and modes derived from a pair of x-ray structures for 20 proteins [7].

3.3 Essential Dynamics and PCA

3.3.1 PCA of Structural Ensembles

A principal component or essential dynamics analysis may be carried out on a molecular dynamics trajectory or any other structural ensemble. It typically consists of three steps. First, the configurations from the ensemble must be superposed, to enable the filtering of internal motions from overall rotation and translation. This is usually accomplished by a least-squares fit of each of the configurations onto a reference structure (see Note 7 ). Second, this “fitted” trajectory is used to construct a variance—covariance matrix that is subsequently diagonalized. The variance—covariance matrix is a symmetric matrix containing, as elements, the co-variances of the atomic displacements relative to their respective averages for each pair of atoms for the off-diagonal elements and the variances of each atom displacements along the diagonal. Atoms that move concertedly give rise to positive covariances, whereas anticorrelated motions give rise to negative entries. Noncorre-lated displacements result in near-zero covariances (see also Note 8 ). Diagonaliza-tion of this covariance matrix yields a set of eigenvectors and eigenvalues, which are usually sorted such that the eigenvalues are in decreasing order. The eigenvalues represent the variance along each of the corresponding collective modes (eigenvectors) and usually a small number of modes suffice to describe the majority of the total fluctuation. As a third step, the original trajectory may be analyzed in terms of the principal components. To this end, the trajectory is projected onto each of the principal modes to yield the time behavior and distribution of each of the principal coordinates (see also Note 9 ). Often, two- or three-dimensional projections along the major principal components are used to allow a representation of the sampled distribution in configuration space or to compare multiple ensembles along the principal modes of collective fluctuation. These projections onto single or multiple principal coordinates can also be readily translated back into Cartesian space to yield an ensemble or animation of the motion along a selection of principal coordinates.

In contrast to standard NMA, a PCA can be carried out on any subset of atoms, and, for proteins, usually only C^α or backbone atoms are taken into account (see also Note 9 ). This has the advantage that the storage and diagonalization of the covariance matrix is less demanding, whereas the main collective modes are very similar to an all-atom analysis [9, 26]. An additional advantage of a backbone-only analysis is that artificial apparent correlations between slow side-chain fluctuations and backbone motions are not picked up by the analysis. A PCA may be compared with results from a standard NMA. However, to this end, one must perform an all-atom PCA and the fluctuations must be calculated from mass-weighted displacements [21]. This form of PCA is often referred to as “quasiharmonic analysis.” If an all-atom analysis is required, an approximation may be used to retrieve only the principal modes of fluctuation, that alleviates the need to store and diagonal-ize the full matrix [26]. As mentioned above, the PCA technique is not limited to the analysis of molecular dynamics trajectories but can be carried out on any ensemble of structures. It can, e.g., be carried out to derive the principal modes from sets of x-ray structures [27], to compare simulation data with experimental conformations [28–30] (see also Fig. 2), or to derive search directions from multiple homologous structures to aid homology modeling [31].

3.3.2 Convergence of PCA Results Derived from Molecular Dynamics Simulations

Principal components derived from different simulations or simulation parts allow us to compare the major directions of configurational space and sampled regions and to judge similarity and convergence. It has been observed that sub-nanosecond protein molecular dynamics simulations suffer from a significant sampling problem, resulting in an apparently poor overlap between the principal components extracted from multiple parts of these trajectories [32, 33]. Nevertheless, it was observed that despite the fact that individual principal components may be different, the subspaces that are spanned by the major principal components converge remarkable rapidly and show a favorable agreement not only between different simulation results, but also between simulation and experiment [28, 30,34,35], see also Fig. 2.

The anharmonic dynamics along the principal modes of collective fluctuation that corresponds to the jumping between multiple local energy minima results in a diffusive dynamics of the principal coordinates [16, 17]. The analogy of this diffusive dynamics to a multidimensional random walk allows one to assess the convergence of the dynamics along the principal (and usually slowest) modes by comparison of the time evolution of the principal coordinates with cosines that would result from random diffusion [36,37]. A high cosine content typically indicates a nonconverged trajectory. Note, however, that a lack of convergence of the dynamics along a set of modes does not necessarily also imply that the directions of such modes or the subspace they span are not converged or poorly defined. Provided that a sufficiently converged trajectory is available, thermodynamic properties may be derived as ensemble averages and can be readily mapped onto the principal coordinates to yield, e.g., free energy landscapes (see Fig. 3).

3.3.4 Enhanced Sampling Techniques

Knowledge of the major coordinates of collective fluctuations opens the way to develop specialized simulation techniques tailored toward an efficient or even systematic sampling along these coordinates, thereby alleviating the sampling problem inherent to virtually all common computer simulations of biomolecular systems today. The first attempts in this direction were aimed at a simulation scheme in which the equations of motion were solely integrated along a selection of primary principal modes, thereby drastically reducing the number of degrees of freedom [9]. However, these attempts proved problematic because of nontrivial couplings between high-and low-amplitude modes, even though, after diagonalization, the modes are linearly independent (orthogonal). Therefore, instead, a series of other techniques has prevailed that takes into account the full-dimensional simulation system and enhance the motion along a selection of principal modes. The most common of these techniques are conformational flooding [11] and essential dynamics sampling [12–14]. In conformational flooding, an additional potential energy term that stimulates the simulated system to explore new regions of phase space is introduced on a selection of principal modes (Fig. 4), whereas, in essential dynamics, sampling a similar goal is achieved by geometrical constraints along a selection of principal modes. More recently, the concept of conformational flooding was reformulated in the context of metadynamics [39]. These techniques have in common that a sampling efficiency enhancement of up to an order of magnitude can be achieved, provided that a reasonable approximation of the principal modes has been obtained from a conventional simulation.