# Atomic Torsional Modal Analysis for high-resolution proteins

###### Abstract

We introduce a formulation for normal mode analyses of globular proteins that significantly improves on an earlier, 1-parameter formulation tirion96 () that characterized the slow modes associated with protein data bank structures. Here we develop that empirical potential function which is minimized at the outset to include two features essential to reproduce the eigenspectra and associated density of states in the – frequency range, not merely the slow modes. First, introduction of preferred dihedral-angle configurations via use of torsional stiffness constants eliminates anomalous dispersion characteristics due to insufficiently bound surface sidechains, and helps fix the spectrum thin tail frequencies ( – ). Second, we take into account the atomic identities and the distance of separation of all pairwise interactions, improving the spectrum distribution in the – range. With these modifications not only does the spectrum reproduce that of full atomic potentials, but we obtain stable, reliable eigenmodes for the slow modes and over a wide range of frequencies.

###### pacs:

87.15.H-, 87.15.ad,Proteins, nanomachines of the living world, provide a bewildering array of inter and intra-cellular services: catalysis, transport, scaffolding, support, cell-signaling and replication among others. Each protein functions with extraordinary almost paradoxical efficiency, specificity and accuracy. A single molecule of acetylcholine esterase, for example, breaks down 25,000 neurotransmitters molecules per second, a rate that seems to be diffusion-limited.

Proteins are linear polymers of 20 different amino acid building blocks whose sequence specifies all covalent bonds. The so-called primary sequence of over 8 million proteins are known, and vary in length from a few dozen to tens of thousands of amino acids. The primary sequences fold into distinct secondary patterns of helices, sheets and loops that in turn combine to create distinct globular arrangements that range in cross-sectional diameter from ten to hundreds of Angstroms. X-ray crystallographers have determined the atomic arrangements of nearly 100,000 protein molecules, revealing recurring folding patterns and structural motifs PDB ().

Proteins are not rigid molecules. Ever since the first crystal structures, of hemoglobin and myoglobin, revealed the absence of entry and exit channels for ligand to reach the buried heme binding site perutz (), motility was understood to be vital to function frauenfelder (). Protein molecules possess motility spectra over time-scales of at least 14 orders of magnitude, from the femtosecond regime for bond breaking, to the pico- and nanosecond regimes for events like flipping of interior and surface indole and benzene rings, to the micro- to millisecond regimes for conformational changes such as occur during allosteric re-arrangements of secondary and tertiary structural elements, to millisecond to seconds for rigid-like tumbling. As each of these dynamic processes is best characterized by different physical theories such as quantum vs. classical mechanics, diffusion and Brownian motion, characterizing and sorting the many dynamic features of proteins is challenging.

Ingenious methods to track local and global motions experimentally continue to be developed using X-rays, NMR spectroscopy, and inelastic and quasi-elastic neutron scattering, to name a few kern (). Experiments detect complex net motility spectra of proteins in particular environments, and isolating specific molecular events from the data remains difficult. Theory segregates events and models catalytic, diffusive and allosteric events separately into fast local and relatively slower nonlocal motions.

Rapid, local catalytic events are well-studied and modeled using quantum mechanical and detailed semi-empirical formulations warshel (). Extended, global motions are less well characterized as they are hard to isolate experimentally and require lengthy simulations. One attractive candidate to characterize global motion is normal mode analysis (NMA) where instantaneous, equilibrium motility spectra are computed go (); karplus (); levitt85 (). The highly idealized condition of an absence of unbalanced forces permits rapid and simple snapshots of the internal flexibilities of proteins, unattainable by standard MD, and provides fascinating insights into protein architecture and design.

Do proteins ever exhibit the motility spectra suggested by NMA? Experimental studies, from Mossbauer spectroscopy, inelastic neutron scattering and X-ray crystallography, demonstrate the existence of equilibrium or harmonic motion in proteins only below a critical “glass transition” temperature of around 180-220K petsko (). Do the eigenmodes then retain relevance at physiological temperatures? At least two lines of evidence suggest they do. High resolution X-ray crystallography provide isotropic and sometimes anisotropic Debye-Waller values for each scattering center. These experimental values, a reflection of the net mobility of each atom in the crystalline environment, are surprisingly well-modeled by the normal mode spectrum of an isolated molecule. Second, and equally surprising perhaps, is the successful use of normal modes to model non-equilibrium, allosteric transitions between known, crystal structures of apo- and holo-forms of the same enzyme petrone (); bray (); omori (); wolynes (). Both lines of inquiry suggest that the very simple and idealized normal mode computation can provide helpful insights above the glass transition temperature as well.

Normal mode analyses provide different levels of detail depending on the formulation adopted in setting up the algorithm. The use of atomic coordinates vs. reduced coordinates, internal Cartesian vs. dihedral coordinates, simple vs. detailed energy potentials, as well as the resolution and stability of the coordinate file all affect the computed eigenfrequency spectra and associated amplitudes as well as the appearance of each independent mode, over the entire frequency range.

Since demonstrating that an idealized, one-parameter energy potential reproduces very well the overall appearance of temperature factors and therefore the slow modes that contribute primarily to this distribution tirion96 (), the question remained: What details were sacrificed? Might it be useful to more faithfully reproduce eigenfrequency spectra observed by spectroscopy and standard NMAs that use full, molecular mechanic potentials and exact molecular topology na (); song ()? Is it indeed possible, within the Hookean approximation of tirion96 (), to achieve the level of detail seen in standard NMA? If yes, we could (a) begin to correlate observed spectral characteristics with specific structural features, and (b) obtain more realistic flexibility descriptions of proteins whose structures are known.

The current work represents an effort to address these questions. We demonstrate that use of (1) complete atomic coordinate files, (2) dihedral coordinates, and (3) more realistic energy formulations than single-parameter Hookean springs can more realistically reproduce experimental and theoretical vibration spectra. Dubbed ATMAN for Atomic Torsional Modal ANalysis, this method provides a rapid visual assessment of the internal flexibility characteristics inherent in the overall design of each protein and may be useful to study design characteristics that contribute to observed maxima in spectroscopic data.

The general formalism for NMA was presented in levitt85 (); go (); karplus (). Briefly, given a protein consisting of atoms of mass , located at , and a potential energy function , one selects a proper set of generalized degrees of freedom and solves the co-diagonalization problem

(1) |

where the elements of the force matrix (also known as the Hessian), and the mass-weighted inertia matrix , are:

(2) |

and the potential is assumed to be minimized at . The eigenfrequencies are given by the elements of the diagonal matrix ; , and the corresponding eigenmodes are the columns of , normalized according to . With this normalization, the amplitude of mode , resulting from thermal activation at temperature , is . Thus, the dominant contribution to the total mean-square fluctuation, given by the sum over all the modes, comes from the low-frequency modes.

The various NMA approaches differ mainly in their choice of the degrees of freedom and the form (and level of detail) of the potential energy function, . Especially for proteins whose structure files are incomplete or known to low resolution, a suitable approach is to select the Cartesian coordinates of reduced masses centered at the carbons, and select an suitable for such a coarse-grained description eyal (); phillips (). For high resolution structures, coarse-graining of the atomic coordinates is unnecessary.

We strongly advocate the use of dihedral degrees of freedom to compute slow modes of proteins. When comparing two protein sequences, one in an active, folded form and one in an unfolded, extended form, very little variation in bond lengths and bond angles exists. The two forms vary mostly in the distribution of 4-atom rotational torsional degrees of freedom (for this reason, algorithms to predict protein folding using dihedral angles abound singh ()). We use the “soft” torsional angles of proteins: all main-chain (MC) and as well as side-chain (SC) angles (except for the MC angles of prolines) as detailed in the L79 potential of Levitt L79 (). Not only does this reduce the size of the computation compared to atomic coordinates, but more importantly use of torsional degrees of freedom necessarily maintains stereochemically correct bond lengths and bond angles throughout.

A suitable potential energy is needed to quantify the energetic costs of distorting the protein’s configuration along each internal coordinate . Detailed semi-empirical formulations such as ENCAD ENCAD () and CHARMM CHARMM () exist that accurately model the various types of relevant interactions: van der Waals, hydrogen bonds, etc. A drawback of these potentials for the current use is that the protein configuration as read from the PDB, , is not initially at a local minimum, as required by a NMA. Minimizing a complex potential function of hundreds of parameters is problematic: The potential landscape consists of many close local minima, making it difficult to compare among the results from different studies; and the minima one finds numerically exhibit various unstable modes, due to inaccuracies, casting doubt on the frequencies and shape of the remaining, stable modes.

In previous work tirion96 () we proposed that a potential centered around the PDB configuration of the protein gives a reasonable description of the modes. Specifically, we suggested

(3) |

where the sum is restricted to non-bonded atom pairs separated by at least three bond lengths and a certain interaction range, (see below). Each pair-wise potential term followed a simple Hookean form:

(4) |

where , and is the initial separation between each interacting -pair as read from the PDB file. This potential has the obvious advantage of being at a well defined (stable) minimum at the outset.

In tirion96 () we also argued that for the slow modes of densely packed globular proteins one can use a single universal value for the spring constants , since they depend on a large number of non-covalent interactions (NBI) whose sum approaches a universal form, governed by the central limit theorem, regardless of the details of individual pair-wise potentials. The actual value of the universal spring constant was adjusted to fit other NMA work. We now discuss how to generalize this single bond-strength formulation, and replace with atom- and distance-dependent non-bonded stiffness constants, , as well as add a new stiffness constant, , for the dihedral angles. The advantages, as we will see, include more realistic frequency spectra and a greater range of applicability due to the greater stability of side chains that may be unbound or minimally bound in the PDB coordinate file.

The general principle is as follows. Given the PDB coordinates of a protein, , and a reliable empirical potential function , we assume that the minimum of lies nearby summa (). In that case, one can approximate the Hessian , from Eq. (2), by simply evaluating the second derivatives at instead of at the minimum of . The fact that first derivatives are not zero at is not a problem because they do not contribute to the Hessian. This can be seen by Taylor-expanding pair-wise potentials around the PDB coordinates: . The contribution of to vanishes if the relation between the coordinates and the generalized degrees of freedom is linear, , as coded in the matrix. Since we are replacing the potential with second-order terms, ignoring first derivatives, this guarantees an extremum point. To ensure that the point is a minimum and the PDB configuration represents a stable equilibrium we exclude the domain where .

Consider, for example, a van der Waals (vdW) force between atoms and given by the Lennard-Jones 6-12 potential

(5) |

where and are the potential well-depths and the distance between and at the minimum of , respectively, while is the distance between the atoms. Note that need not necessarily equal the starting value from the PDB file, . For ’s contribution to the Hessian, it suffices to consider only the second-order derivative:

(6) |

This expression turns negative for , and we must exclude that domain in order to avoid potential instabilities. Chopping off the domain results, however, in an effectively reduced interaction range between the atoms. In practice, empirical potentials consider interactions between atoms and up to a certain range that is computed by adding an arbitrary cutoff distance to the sum (or a multiple of the sum) of the van der Waals radii:

(7) |

and we would like to extend our interaction range all the way to as well. This can be done by “stretching” our domain. We first note that admits a simpler (but approximate) form:

(8) |

see Fig 1.

The stretching of the interaction range is then achieved by mapping , to yield

(9) |

This atom- and distance-dependent spring strength replaces the universal constant of Eq. (4).

As a second example consider the potential function for the rotation of the dihedral angles coordinates. Typically:

(10) |

where is the potential depth, is the number of minima along the rotation coordinate , and is the value of at the potential’s minimum. Once again, if we are using dihedral angles as our internal coordinates, only the second derivative contributes (in this case, to the entry):

(11) |

To ensure stability, we limit ourselves to the domain where . The reduced interaction range that results from this procedure can be stretched from to a larger arbitrary cutoff :

(12) |

Each of the torsional degrees of freedom can then be endowed with a potential , where is computed as above, with the appropriate values of , and prescribed by the empirical parent potential. In practice, for our purposes of fitting the total distribution of modes to that of a detailed potential, we find that it is sufficient to assume a universal constant for all torsional degrees of freedom (instead of the more detailed expression (12)) and to take the range of to include the whole circle.

It should be clear from these examples that any potential formulation, whatever its complexity and degree of sophistication, can be treated similarly. We note that our procedure of cutting off the range of the variables when 2nd derivatives become negative (unstable) is not actually necessary in a properly minimized configuration. Without this procedure, our analysis leads to instabilities, and the cutting of the range (and extending and stretching to some additional arbitrary distance, to compensate) is the price one pays for approximating the real minimum by the PDB configuration.

While our formulation of the potential energy is still oversimplified, it retains the advantage of defining the initial configuration to be an exact energy minimum, obviating the need for structure-distorting energy minimizations and the concomitant risk of negative eigenvalues due to incomplete minimizations. The simple potential of tirion96 () yielded slow modes (below the range) that were comparable to those found using more sophisticated, detailed potentials. However, for the distribution of the modes deviated markedly from that obtained by standard methods (see curve (f) in the inset of Fig. 3). In the following we show that making the single, non-bonded stiffness strength from tirion96 () atom-type and distance dependent, as in Eq. (9), plus adding a simple dihedral potential with the same constant for the various torsional degrees of freedom, achieves a better description of the distribution of modes throughout the whole range of frequencies (up to ). Interestingly, a recent independent study, focusing on B-factors rather than on the distribution of modes, by H. Na and G. Song na () reached a very similar conclusion: Including detailed van der Waal interactions and a universal dihedral potential achieves the best fit to standard NMA in their excellent and comprehensive case study.

To demonstrate, we computed the normal modes of xylanase 10A (PDB entry 1GOK confomer A), a 302 amino-acid enzyme with a mass of nearly 33,000 daltons produced by a thermophilic ascomycetous fungus. Xylanase 10A, a member of the family 10 glycoside hydrolases, cleaves internal bonds of xylan, a major polysaccharide component of plant cell walls, releasing shorter xylose subunits. The crystal structure of 1GOK has been determined to resolution leggio () and is seen to possess a classic TIM barrel fold, Fig. 2, consisting of eight parallel -strands forming a central tube connected by 8 -helices external to the -strand tube. The structure has been likened to a ‘salad bowl’ as one face has a larger radius, of about , while the opposite face has a smaller radius of about collins (). The active site consists of a substrate-binding cleft running the length of the larger face of the ‘salad bowl’ and includes two conserved residues, Glu-131 and Glu-237, that act as the catalytic proton donor and nucleophile for the hydrolysis. The crystal structure reveals that three conserved tryptophans, at locations 87, 267 and 275, form an aromatic cage around the active site.

To compute the normal modes associated with this structure, we built in MC amide hydrogens reduce () to obtain a 2598-atom system with 589 MC and and 471 SC dihedrals, resulting in a total of 1060 degrees of freedom. We computed the F and H matrices in Eq. (2) using both analytic and numeric algorithms to check for consistency. The use of non-Cartesian degrees of freedom necessitates elimination of overall translations and rotations about the center of mass in the computation of . Failure to do so results in incorrect frequencies and eigenmodes, including zero- or numerically unstable modes.

Our original formulation of the Hookean NMA tirion96 () set the dihedral stiffness constants to zero. Compared to the cumulative non-bonded interactions of thousands of atom pairs, the contribution to the total energy of each dihedral was considered negligible. Careful examination of the eigenspectra, however, revealed that the slowest modes often included ones resulting in root mean square (RMS) displacements discontinuous with neighboring modes. Studying the animations of these anomalous eigenmodes revealed isolated surface sidechains sweeping wide arcs, independent of the remainder of the molecule. These side chains were characterized by relatively small numbers of interatomic interactions to stabilize their conformations, which resulted in the anomalous RMS displacements. In past studies, we replaced such residues with alanine residues lacking sidechain degrees of freedom, obviating the problem. However, introducing physically plausible constraints on the dihedral bond deformations removes the need to alter the PDB coordinates. We found that a value of stabilized all surface side chains and resulted in reliable and accurate eigenspectra.

For the non-bonded, vdW forces we used Eq. (9) with , to identify all non-bonded atom pairs separated by at least 3 bond lengths and within an interaction distance given by Eq. (7). With values ranging from 1.5 to 2.4 , this resulted in 8,000 to 21,000 non-bonded interactions. The resulting stiffness constants were computed using Eq. (9) and set to . is a constant introduced to adjust the overall scale of the eigenfrequencies and was chosen, in each case, to make the frequency of mode 1 equal to . This required values ranging from 160. to 4.0. The values of the van der Waals radii, , and the Lennard-Jones well depths were taken from ENCAD ENCAD ().

With the various van der Waals radii and determined from standard potentials, there remain only three parameters to adjust: The cutoff distance ; the overall constant , in the relation ; and the dihedral curvature constant . The distribution of modes of globular proteins (at all frequencies), as obtained from standard potentials, tend all to fall under one universal curve tirion93 (); universality () and we use this curve to fit our parameters. The universal curve provides two important clues, as and and have very different effects on the shape of the distribution. This, coupled with the fact that the slowest mode of such a protein as xylanase is expected to be around constrains the value of the three parameters quite effectively remark ().

The resulting eigenfrequency spectra, , are presented in Fig. 3 as histograms with bin size. The spectra (a-d) for different numbers of NBI are superposed on the eigenvalue spectra of various proteins computed by Levitt levitt85 () using a standard energy formulation and energy-minimized coordinates. In an effort to better characterize the effects of the current formulation, we include in the inset of Fig. 3 equivalent spectra of 1GOK using the earlier formulation presented in tirion96 ()(f), as well as that formulation combined with the constraints on the dihedrals (g). In all cases, the low-frequency modes overlap, indicating that the various formulations identify the character of the slowest modes equally well. The simplest Hookean formulation (f), using identical spring-constants for all non-bonded interactions and no dihedral constraints, obtains a sharp peak and narrow distribution centered around , with no high-frequency tail. Interestingly, inelastic neutron scattering experiments by Cusack and Doster identify such a well-defined maximum in the dynamic structure factors at around in myoglobin for temperatures below 180K cusack (). This peak is well resolved in all NMA formulations we tested.

Inclusion of dihedral stiffness constraints adds a second maximum to these spectra, as seen in (g), around , widening the distribution and introducing higher-frequency, tail-end contributions. While higher frequency contributions are observed, the additional peak at is not reported on in experiments and not resolved by standard NMA, and indicates that the formulation remains overly simplistic. Use of scaled non-bonded stiffness constants adjusted the spectra significantly (a-d). Decreasing the number of NBI has the effect of widening the distribution over a larger frequency range and de-emphasizing the peak. The largest number of NBI tested, at 21,000 [curve (a)], produces a distribution similar to that of (g) that used a single stiffness constant, for all interacting atom pairs. At equal to 1.6Å, resulting in 9000 nonbonded interactions, we obtained the best fit to the theoretical curve of Levitt. Decreasing the cutoff further tends to spread the distribution even more evenly, as seen in (d) for .

Studying each of these cases (a-e) via animations as well as RMS fluctuations demonstrates an interesting pattern, a pattern further clarified in the limit of zero non-bonded interactions. In the case of zero NBI and only dihedral stiffness constraints, we obtain an eigenspectrum with purely one-dimensional traits: a broad distribution over the entire frequency range, and a density of states for small frequencies. As the number of NBI increases, the spectrum begins to shift toward that of a two-dimensional solid, with distinct peaks and with . In particular, for cutoff values that permit fewer than about 5,000 non-bonded interactions, the spectra remain largely two-dimensional and no “hitching” of neighboring (non-bonded) regions toward cohesive motion is observed. As the number of NBI increase, neighboring strands begin to “hitch” together and move cohesively. Increasing the NBI still further strengthens the cohesion and larger domains move in lockstep. The largest cutoff values lead to a motility spectrum exhibiting very smooth motions more typical of three-dimensional solids with .

Larger cutoffs resulting in more cohesive motions appeal to the eye in animations due to their ultra-smooth, rubber-like deformations. However, the computed frequency spectra indicate that this regime does not best model either spectroscopic or theoretical results, which indicate a broader maximum and softer molecule, even in the case of the temperature regime below 180K. We find that small cutoff distances leading to minimal cohesion, or “hitching,” best reproduce eigenspectra characteristics. This condition also ensures sensible activation amplitudes relevant to computing mean square displacement values and describe stereochemically sound structural deformations. (What constitutes “minimal cohesion” seems to depend on the resolution of structures, with lower-resolution structures obtaining larger cutoff values than high-resolution structures.)

Animated GIF sequences of the first three slowest modes are available at gifs (). These sequences were created by PyMol and rendered in GIF format by LICEcap, and use a surface rendering mode to simplify the representation of the many atoms. In each left-hand panel, the structure is oriented as in Fig. 2 so as to look at the top of the “salad bowl” where the ligand-binding groove and active site are situated. The active site residues, Glu-131 and Glu-237, are colored magenta, while the active-site cage formed by the triad Trp-87, Trp-267 andTrp-275, is colored green. The ligand-binding groove runs horizontally, with subsites binding nonhydrolyzed subunits binding to the left of the active site and subsites holding hydrolyzed reactants to the right. In mode 1, the molecule was rotated around a vertical axis by to create the right-hand panel. In mode 2, the molecule was rotated around a vertical axis to create the view in the right-hand panel looking down from the top of the molecule to demonstrate the twisting nature of the motion. In mode 3, the molecule was rotated around a vertical axis to see the bottom of the molecule in the right-hand panel.

What is immediately apparent is that while the computation of eigenmodes proceeds without any information about ligand-binding propensities and active site residues information, the slowest 3 modes each center on precisely these regions. It would seem that structure has been selected to promote dynamic access to the active site and ligand-binding groove. The slowest mode pertains to an up-and-down chewing type motion, the second slowest mode resembles a rotational grinding motion, while the third slowest mode describes a side-wise rolling type motion. These eigenmodes appear ideally suited to enhance entry of ligand and release of products from the binding groove, and to possibly influence the translocation of the polysaccharide chain along the binding-groove. While detailed, mechanistic models may be unwarranted by these simple calculations and visualizations, it would appear that protein architecture co-evolved with catalytic functionality, precisely in order to exploit innate flexibility characteristics to optimize yield.

To better appreciate the combination and distribution of dihedral angles contributing to each mode and their associated domain motions, we overwrote the Debye-Waller temperature factor, , in the PDB coordinate file with normalized and binned dihedral values for each eigenmode vector. To visualize the distribution of the dihedral strengths for each mode, we associated all the torsions with their associated MC amide nitrogens and all the torsions to the MC carbonyl carbons.

Studying the distribution of dihedral shifts for mode 1, for example, explains how the C-terminal and N-terminal domains move en-masse, and how the overall appearance is of two domains that come together across the ligand-binding groove, as in the first animation gifs (). angles in mode 1 with normalized values of 0.5 or higher include Gly-41, Pro-45 and Pro-200. angles with values over 0.5 include Gly-41, Ile-201 and Thr-208. These “hinge-points” permit the relative reorientation of regions, with the N and C terminal regions moving as a unit above the binding-groove. Gly-41 at the N-terminal end of -strand 8 is seen to move in phase with the upper domain, while Pro-45 at the C-terminal end of the same strand moves in phase with the lower domain. Likewise, Pro-200 preceding -strand 4 moves in phase with the lower domain, while Thr-208 at the C-terminal end of the same strand moves along with the upper domain. So mode 1 would seem to have an upper domain extending across residues 1-45 as well as residues 208-302, while the lower domain seems to extend from residues 46-200.

In this manner we tried to isolate “hot spots” that enable this motion and sought to mutate these hinge-points to suppress the observed motility pattern. Suspecting Gly-41, Pro-45, Pro-200 and Thr-208 as essential to the motion inherent in mode 1, we “mutated” these residues by eliminating their associated MC dihedral degrees of freedom (reducing the size of the computation by 8 degrees of freedom, in other words). Suspecting we might suppress mode 1 in this way, we wondered if mode 2 would now make the dominant contribution to the RMS fluctuations. In fact, this did not happen. Mode 1 reappeared as before, but with a different distribution of hinges, to once again present as two domains that come together across the binding-groove. And this seems to be a general finding: the slow modes are a characteristic of the overall architecture of the molecule, independent of any particular residues (publication in preparation).

In conclusion, we find that NMA does support efforts to characterize global protein motion on relatively slow time scales. A straightforward computation that characterizes the full motility spectrum at a single instant, NMAs require neither excessive memory nor CPU capabilities (current computations were run on a Mac Mini with a single Intel Core i5 processor). Visualizations of eigenmodes via animations, using the many tools available for rendering atomic coordinate files, are effortless. A simple sequence of 5 to 9 PDB entries, flashed consecutively on the screen, creates a dynamic effect with a level of realism not attainable using reduced coordinates.

Furthermore, features of the eigenfrequency spectra contain additional clues to structural characteristics. Cusack and Doster, for example, deduced that the experimentally observed maximum at around is not due to internal deformations of secondary structures, but rather to cross-chain interactions, such as the relative motions of helices cusack (). Our current analyses support this interpretation, as a reduction in the non-bonded cutoff range eliminates the “knitting together” of non-bonded regions and the pronounced peak at . Furthermore, we find the appearance of a second peak, near , by the inclusion of dihedral stiffness constants. The height and width of this additional peak depends sensitively on the value(s) of chosen, as does the appearance of the resultant high-frequency tail-end of the distribution. Here we used a small value, equal to across all dihedrals, simply in order to stabilize the structure. However, it is certainly possible to adopt more realistic values for each dihedral angle in order to more faithfully reproduce the spectra.

Inclusion of dihedral stiffness constants stabilizes surface side-chain residues and C- and N-terminal residues that typically obtain less packing constraints than inner residues, and enable more reliable diagonalizations, ones that do not obtain anomalous slow frequencies. However, these spectra still lack the width typical of standard NMA. In order to reproduce this characteristic, we needed to soften or lessen the energetic contribution of nonbonded atom pairs at greater distances. Cutting the nonbonded cutoff distance to the smallest permissible while still maintaining minimal cohesion between nonbonded strands is helpful in softening the distribution but does not appear to be sufficient. We find that de-emphasizing the contributions of pairwise interactions at larger distances recreates the soft, wide distribution successfully. Dividing the total numbers of NBI by the number of atoms in the molecule provides a crude estimate of the numbers of nonbonded interaction per atom, on average. Values over 12/atom lead to stiff molecules with narrow distributions and steeper slopes at small frequencies. Optimal results occur when the average numbers of interactions per atom are no more than 6, which would seem to include only nearest neighbor interactions, within a single “shell”. Interestingly, at these smaller cutoff values, a shift in the distribution of types of nonbonded interactions occur. If we divide the list of non-bonded interactions into 3 groups, for those that pertain to MC-MC, MC-SC and SC-SC interactions, the largest group depends on the cutoff value. In particular, at small cutoff values, MC-MC type interactions predominate, while as the cutoff increases, MC-SC type interactions become dominant, perhaps providing further rationalization of how linear chains folded into compact shapes obtain spectra with spectral dimension greater than 1 but less than 2 elber ().

How might solvent affect computed modes? Frequency shifts due to solvent depend on whether proteins are overdamped or underdamped in their aqueous environments. Oscillators in underdamped environments experience no shifts in their eigenfrequencies while overdamped oscillators do. Experimental evidence favors the underdamped situation for proteins. Recently, for example, Turton and collaborators reported results with extremely sensitive optical Kerr-effect spectroscopy of lysozyme turton (). They found global vibrational modes that are underdamped, and concluded that “the ligand-binding coordinate in proteins is underdamped and not simply solvent-controlled as previously assumed” turton (). The current analysis adds to this intriguing finding by demonstrating motility patterns centered on the active site.

We conclude that slow (ATMAN) modes of PDB entries can provide novel insights to crystallographers and biochemists who aim to elucidate structure and function correlates. While “a picture is worth a thousand words,” we suspect an animation is worth another hundred or so and foresee a future where each PDB entry may obtain another tab for visualizing the consequence of its particular architecture by gifs demonstrating the first few slowest modes, and with the possibility of downloading such sequence files to study that information independently wako (); song ().

###### Acknowledgements.

MMT gratefully acknowledges support from Michael Levitt who emailed every two or three years to assure us that (a) modes were still very much in vogue and (b) were really easy to do nowadays. Without his continued support my interest would not have been reawakened, and I cannot sufficiently express my gratitude for this. MMT also thanks M. Hildred Blewett whose generous bequest provided the means to re-activate this research. This work is supported in part by the M. Hildred Blewett Fellowship of the American Physical Society, http://www.aps.org.## References

- (1) F. C. Bernstein, T. F. Koetzle, G. J. Williams, E. E. Meyer Jr., M. D. Brice, J. R. Rodgers, O. Kennard, T. Shimanouchi and M. Tasumi, “The Protein Data Bank: A computer-based archival file for macromolecular structures,” J. Mol. Biol. , 112 535 (1977).
- (2) M. F. Perutz and F. S. Mathews, “An x-ray study of azide methaemoglobin,” J. Mol. Biol. 21, 199-202 (1966).
- (3) H. Frauenfelder, S. G. Sligar, and P. G. Wolynes, “The energy landscapes and motions of proteins,” Science 254, 1598-1603 (1991).
- (4) K. Henzler-Wildman and D. Kern, “Dynamic personalities of proteins,” Nature 450, 964-972 (2007).
- (5) A. Warshel, “Multiscale modeling of biological functions: From enzymes to molecular machines,” Angewandte Chemie 53, 10020-10031 (2014).
- (6) N. Go, T. Noguti and T. Nishikawa, “Dynamics of a small globular protein in terms of low-frequency vibrational modes,” PNAS 80, 3696-3700 (1983).
- (7) B. Brooks and M. Karplus, “Harmonic dynamics of proteins: Normal modes and fluctuations in bovine pancreatic trypsin inhibitor ,” PNAS 80, 6571-6574 (1983).
- (8) M. Levitt, C. Sander and P. S. Stern, “Protein normal-mode dynamics: Trypsin inhibitor, crambin, ribonuclease and lysozyme,” J. Mol. Biol. 181, 423-447 (1985).
- (9) D. Ringe and G. A. Petsko, “ ‘The glass transision’ in protein dynamics: What it is, why it occurs, and how to exploit it,” Biophys. Chem. 105, 667-680 (2003).
- (10) J. K. Bray, D. R. Weiss and M. Levitt, “Optimized torsion-angle normal modes reproduce conformational changes more accurately than Cartesian modes,” Biophys. J. 101, 2966-2969 (2011).
- (11) P. Petrone and V. S. Pande, “Can conformational change be described by only a few normal modes?,” Biophys. J. 90, 1583-1593 (2006).
- (12) S. Omori, S. Fuchigami, M. Ikeguchi and A. Kidera, “Linear response theory in dihedral angle space for protein structural change upon ligand binding,” J. Comp. Chem. 30,2602-2608 (2009).
- (13) O. Miyahita, J. N. Onuchic and P. G. Wolynes, “Nonlinear elasticity, proteinquakes, and the energy landscapes of functional transitions in proteins,” PNAS 100, 12570-12575 (2003).
- (14) M. M. Tirion, “Large amplitude elastic motions in proteins from a single-parameter, atomic analysis,” Phys. Rev. Lett. 77, 1905 (1996).
- (15) H. Na and G. Song, “Bridging between normal mode analysis and elastic network models,” Proteins 82, 2157 (2014).
- (16) H. Na and G. Song, “Conventional NMA as a better standard for evaluating elastic network models,” Proteins 83, 259-267 (2014).
- (17) E. Eyal, L. W. Yang,and I. Bahar, “Anisotropic Network Model systematic evaluation and a new web interface,” Bioinformatics 22, 2619-2627 (2006).
- (18) D. A. Kondrashov, A. W. Van Wynsberghe, R. M. Bannen, Q. Cui and G. N. Phillips, Jr., “Protein structural variation in computational models and crystallographic data,” Structure 15, 169-177 (2007).
- (19) H. Singh, S. Singh and G. P. S.. Raghava, “Evaluation of protein dihedral angle prediction methods,” PLOS One 9 (2014).
- (20) M. Levitt, “Molecular dynamics of native protein. I. Computer simulation of trajectories,” J. Mol. Biol. 168, 595-620 (1983).
- (21) M. Levitt, M. Hirshberg, R. Sharon and V. Daggett, “Potential energy function and parameters for simulations of the molecular dynamics of proteins and nucleic acids in solution,” Comp. Phys. Communication 91, 215-231 (1995).
- (22) B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan and M. Karplus, “CHARMM: A program for macromolecular energy, minimization, and dynamics calculations,” J. Comp. Chemistry 4, 187-217 (1983).
- (23) C. R. Summa and M. Levitt, “Near-native structure refinement using in vacuo energy minimization,” PNAS 104, 3177 (2007).
- (24) L. L. Leggio, S. Kalogiannis, K. Eckert, S. C. M. Teixeira, M. K. Bhat, C. Andrei, R. W. Pickersgill and S. Laresen, “Substrate specificity and subsite mobility in T. aurantiacus xylanase 10A,” FEBS Letters. 509, 303-308 (2001).
- (25) T. Collins, C. Gerday and G. Feller, “Xylanases, xylanase families and extremophilic xylanases,” FEMS Microbiol. Rev. 29, 3-23 (2004).
- (26) J. M.. Word, S. C. Lovell, J. S. Richardson and D. C. Richardson, “Asparagine and glutamine: Using hydrogen atom contacts in the choice of sidechain amide orientation,” J. Mol. Biol. 285, 1735-1747 (1999).
- (27) http://adweb.clarkson.edu/~mmtirion/.
- (28) M. M. Tirion and D. ben-Avraham, “Normal mode analysis of G-actin,” J. Mol. Biol. 230, 186 (1993).
- (29) D. ben-Avraham, “Vibrational normal-mode spectrum of globular proteins,” Phys. Rev. B 47, 14559 (1993).
- (30) This last requirement — the frequency of the slowest mode — can be replaced by a more robust statistic, e.g., the fact that about 14% of the modes fall below .
- (31) S. Cusack and W. Doster, “Temperature dependence of the low frequency dynamics of myoglobin. Measurement of the vibrational frequency distribution by inelastic neutron scattering,” Biophys. J. 58, 243-251 (1990).
- (32) R. Elber and M. Karplus, “Multiple conformational states of proteins: a molecular dynamics analysis of myoglobin,” Science 235, 318-321 (1987). end
- (33) D. A. Turton, H. M. Senn, T. Harwood, A. J. Lapthorn, E. M. Ellis and K. Wynne, “Terahertz underdamped vibrational motion governs protein-ligand binding in solution,” Nature Communication 5, 1-6 (2014).
- (34) H. Wako, M. Kato and S. Endo, “ProMode: a database of normal mode analyses on protein molecules with a full-atom model,” Bioinformatics 20, 2035-2043 (2004).