Evaluation of Protein Elastic Network Models Based on an Analysis of Collective Motions

: Elastic network models (ENMs) are valuable tools for investigating collective motions of proteins, and a rich variety of simple models have been proposed over the past decade. A good representation of the collective motions requires a good approximation of the covariances between the ﬂuctuations of the individual atoms. Nevertheless, most studies have validated such models only by the magnitudes of the single-atom ﬂuctuations they predict. In the present study, we have quantiﬁed the agreement between the covariancestructure predicted by molecular dynamics (MD) simulations and those predicted by a representative selection of proposed coarse-grained ENMs. We then contrast this approach with the comparison to MD-predicted atomic ﬂuctuations and comparison to crystallographic B-factors. While all the ENMs yield approximations to the MD-predicted covariance structure, we report large and consistent diﬀerences between proposed models. We also ﬁnd that the ability of the ENMs to predict atomic ﬂuctuations is correlated with their ability to capture the covariance structure. In contrast, we ﬁnd that the models that agree best with B-factors model collective motions less reliably and recommend against using B-factors as a benchmark.


INTRODUCTION
Elastic network models (ENMs) provide a greatly simplified perspective on collective protein motions, and ever since Tirion's pioneering work 1 a variety of different models have been suggested.The success of coarse-grained ENMs has to some extent been surprising, and their simplicity has provided novel insights into diverse protein mechanisms.Among many applications we find: characterization of the mechanics of large allosteric transitions and of motor proteins, 2−5 studies of Brownian dynamics, 6 characterization of structural response in evolution, 7−10 and applications in structural alignment (e.g, Zen et al.). 11The simplicity of their parametrization and the relatively low computational cost of ENMs make them convenient compared to traditional molecular mechanics approaches and often enable analyses that would otherwise not be feasible due to constraints on computational resources.ENMs have also been adopted for the study of fundamental physical properties of proteins (e.g., Reuveni et al.) 12,13 for which it is not clear that more detailed models offer any advantages.
Tirion computed thermodynamical statistics for proteins using an intrinsic dynamics model with a greatly simplified harmonic potential. 1The potential was constructed to have its minimum at the equilibrium coordinates reported from experimentally determined structures.This work demonstrated the feasibility of predicting collective motions from ENMs, and how computationally expensive structure minimization can be avoided.In the subsequent decade, several researchers incorporated this finding into ENMs developed for normalmode analysis, adding further simplifications such as coarsegraining at the residue level. 14,15Later, refinements to these initial coarse-grained models were suggested, e.g., approximate side-chain models 16 and refined backbone modeling. 17,18implifications in the degrees of freedom have also been suggested.Bahar et al. 9 proposed an isotropic model employing a potential which restrains angular deviations of contacts. 19ately, coarse-grained models in torsional angle space have been proposed. 20,21These models have completely rigid bonds between C α atoms adjacent in sequence, thereby reducing the dimensionality by ignoring these degrees of freedom and essentially preserving the form of the potential.
−25 Comparisons to crystallographic B-factors and conformational variability in NMR ensembles also serve as rough validations of the atomic fluctuations predicted by ENMs.−30 Computationally, validation has been done by comparison with molecular dynamics (MD) simulations of more detailed, anharmonic potentials that have been parametrized semiempirically. 16,31,32To guide practitioners in choosing between models and parametrizations and aid in further model development, we find it important to benchmark suggested ENMs and evaluate the effect of the suggested model refinements.Some effort has already been made in this regard.For example, Micheletti et al. 16 and Rueda et al. 31 each compare a small selection of ENMs against MD simulations, while Riccardi et al. 33,34 treat a selection of models focusing on crystalline structures.Recently, Romo and Grossfield 35 and Leioatts et al. 36 have compared different ENM formalisms to each other, using MD simulations as a benchmarking standard.
In the present work we aim to do a comprehensive comparison of proposed parametrizations for different ENMs applied to isolated (noncrystalline) structures.We prefer benchmarking against more detailed computational models for two reasons.First, ENMs are motivated by the statistical properties of interaction averages in collective motions, 1 and we find it important to use methods that can compare the inter-residue positional covariances describing this collectivity, rather than just single-atom variances.Second, B-factors are dominated by contributions from crystal defects rather than thermal fluctuations, 34,37 and the variability in NMR ensembles is due to incomplete experimental data as well as thermodynamic effects.This makes both these experimental sources for atomic fluctuation data too unreliable for benchmarking.
In the following, we compare a selection of six types of ENMs with fixed parameters by calculating the agreement of MDpredicted and ENM-predicted covariance matrices.To this end we use MD trajectories of seven solvated proteins of different size and oligomerization.We then contrast our approach with benchmarks comparing B-factors and ENM-predicted atomic fluctuations.To disentangle the quality of the benchmark from the sensitivity of the measure used for comparison, we also benchmark ENMs by comparing with MD-predicted atomic fluctuations, using the same measure of similarity that we use for comparing with B-factors.The results obtained allow us to comment on the different relevance of the proposed ENMs for protein modeling and identify some factors that are important to consider for these kinds of models.We also provide results that demonstrate shortcomings of the common approach to validation against B-factors.

MATERIALS AND METHODS
Crystal Structures.To compare predictions of covariance structures from ENMs and MD data, we selected seven structures of varying size, secondary structure and domain compositions, and oligomerization.The investigated structures are tabulated in Table 1.For the GroEL subunit, chain A of the listed structure was used, although this crystal is obtained for the entire 14-mer complex.
Elastic Network Models.The ENMs considered in this study model the protein as a network of nodes located at the positions of representative atoms of each residue, typically one node for each C α position.In all the models these nodes are interacting through a harmonic pair potential that is minimal when the nodes are both at their equilibrium positions.The potential used in all but one ENM can be generalized as where ∥∥denote the Euclidean norm, r is a conformation of the ENM, r 0 is the equilibrium conformation, and r i and r i 0 are threedimensional column vectors describing the position of node i in the respective conformations.The force constant k ij is a parameter controlling the strength of the interaction and is specified differently in the different models.In this work, the equilibrium conformation used to calculate the ENMs is in all cases taken to be an X-ray structure model obtained from the Protein Data Bank (PDB). 45,46dding the potential energy terms for each pair of nodes and expanding the result in a Taylor series around r 0 , the potential for an ENM is expressed as where H is the Hessian of the system and T denotes the transpose.H has the form: As only internal interactions are considered, rotations and translations have no energetic cost, and since moreover the Hessian is evaluated at an energy minimum where the energy gradients are zero, the rank of H is at most 3N − 6, with N the number of nodes in the network.The Boltzmann distribution for a system governed by a harmonic potential is a Gaussian distribution with covariance matrix proportional to the inverse of the Hessian.In the case of the potential specified in eq 2, whose Hessian is singular, the covariance matrix of fluctuations can be calculated as the Moore−Penrose pseudo-inverse of H, denoted with superscript +.We subsequently normalize the covariance matrix by multiplying with the reciprocal of its trace: where Σ is the normalized covariance matrix and tr denotes the trace of its argument.This normalization removes a global amplitude factor from the fluctuations and covariances, which none of the ENMs describes realistically anyway. 37The relative magnitude of atomic fluctuations is predicted by the positional variances of the nodes: where σ i 2 denotes the fluctuation of node i and Σ ii denotes the 3 × 3 submatrix of Σ, containing the positional covariance of the three Cartesian coordinates for node i.The Gaussian network model (GNM) 19 and a recent refinement of it 30 use a different pair potential: The Cartesian components of displacements contribute independently to this potential, so the potential energy terms can be summed separately for each coordinate axis: where R n and R n 0 are vectors of size N describing the position of each node along each of the coordinate axes, and Γ is a N × N matrix with Rearranging eq 6 reveals that the potential is dependent on the angle between the displacements of connected residues and hence not generally zero for rigid body rotations.For translations the potential is zero, so the rank of Γ is at most N − 1. Summing the covariances along each coordinate axis gives the expected values of the inner products of the displacements, which are proportional to the elements of Γ + .The normalized matrix of the expected values of the inner products Z is whose elements are commonly treated as scalar valued covariances.To contrast this to the other ENMs, note that under their respective potentials H ij = ⟨(r i − r i 0 )(r j − r j 0 ) T ⟩, while Z ij = ⟨(r i − r i 0 ) T (r j − r j 0 )⟩.⟨⟩denotes the expected value.The relative magnitudes of atomic fluctuations ζ i 2 is: which has the same interpretation as σ i 2 for the other ENMs.ENM Parameterizations.In the following the six ENMs we have treated will be presented.Apart from the two forms of the potential energy function (eqs 2 and 6), the ENMs are characterized by how the values for the force constants are determined.We have focused on anisotropic C α models in Cartesian coordinates, in addition to the GNM which was included since it is a very popular model.
Uniform Force Constants (ANM and GNM).The two simplest models, the GNM 19 and the anisotropic network model (ANM), 15 define the force constants as (11)   where d ij 0 = ∥r i 0 − r j 0 ∥, γ is a uniform force constant, and c defines a cutoff distance within which nodes are defined to be interacting.We only treat normalized covariance matrices and atomic fluctuations, which are not affected by the parameter γ, so its value is not important.The parameter c of the GNM is generally chosen to include interactions within the first coordination shell around interior residues, often with a value around 0.7 nm. 19For the ANM there is no clear consensus on a transferable value for this parameter.Cutoffs in the range of 0.7 nm can in this case lead to under-restrained degrees of freedom, causing local overestimates of atomic fluctuations. 15In the case of completely unrestrained degrees of freedom implied in protein deformations, this model has a Hessian of rank lower than 3N − 6.General recommendations for the cutoff in the case of ANM vary from the range of 0.8−0.3nm 47 to the range of 1.5−2.4nm. 27For both models it has been suggested that the appropriate cutoff could be chosen by fitting to B-factors. 27,48We will refer to these models subscripted with the cutoff parametrization in nanometers, e.g., ANM 0.8 for an ANM model with c = 0.8 nm.When investigating the ANM we have represented it with one parametrization from either range (ANM 0.8 and ANM 1.8 ), illustrated in Figure 1.
Distance-Dependent Force Constants (HCA and pfANM).The model introduced by Hinsen et al., 6 which we will refer to as the harmonic C α potential model (HCA), determines the force constants as where the parameter c is set to 0.4 nm to separate interactions between C α atoms adjacent in sequence from the other interactions.Long range interactions are proportional to an inverse power of six of the equilibrium distance between the interacting nodes.The parameters a 1 , a 2 , and b were chosen by fitting to an all-atom normal mode calculation for crambin. 6This is also the parametrization we have used here.The parameter-free anisotropic network model (pfANM) 49 defines force constants by the inverse-square of the equilibrium distance between the interacting nodes: On the request of a reviewer we also tested the dependence of pfANM on the choice of exponent, and in particular to the alternative choice of −6 which is very similar to the HCA model.Results are presented in the Supporting Information.Backbone Refinement (REACH).In a series of publications, Smith and co-workers investigated several aspects of coarsegrained modeling of protein interactions with the realistic extension algorithm via covariance Hessian (REACH) model. 17,50,51In particular, refined force constants for backbone interactions were considered.Long range interactions were treated with a fast decaying distance-dependent function.We where s is the sequential distance between the interacting residues, so that s = 1 for adjacent residues.More generally, two interacting residues are separated by s − 1 other residues along the protein sequence.We have used the parameters obtained from fitting to MD simulations of a myoglobin dimer. 17or a visual comparison of the distance dependence of the force constants in HCA, pfANM, and REACH, see Figure S3.
Implicit Side-Chain Representation (βGM).The β Gaussian Model (βGM) 16 uses a force constant definition very similar to that of the ANM and GNM, but sets itself apart from these by modeling side-chains.Nodes representing C α atoms are located at their equilibrium position, and side-chains are represented for all residues except glycines by a node we will refer to as C β .The position of C β is given by where r βi 0 and r αi 0 denote the position of the C β and C α atoms of i, respectively, and η is chosen to be 0.
with x αi and x βi denoting the displacements of the C α and C β atom of residue i, respectively.Substituting eqs 15 and 16 into eq 2 gives a Hessian of C α coordinates incorporating the effects of interactions with C β atoms: H H HD HD D HD () TT (17)   where H αα , H αβ , and H ββ are the Hessians for interactions between C α atoms, between C α and C β atoms, and between C β atoms, respectively.D is a 3M × 3N-matrix encoding eq 16, the transformation from N C α displacements to M C β displacements.In principle, four force constants are defined: b 00 00 (18)   where t(i) denotes the type of the node (C α or C β ).Following Carnevale et al. 8 we set κ αα = κ αβ = κ ββ , κ β =2κ αα and c = 0.75 nm.
Atomistic Sampling Model.From MD trajectories the covariance matrix V is obtained as a sample statistic: where X j is a column vector of dimension 3N specifying the jth conformation in the trajectory, roto-translated to minimize the root-mean-squared distance (RMSD) over all C α atoms to the conformation reported in the PDB.M is the number of conformations in the trajectory, and X ̂is the mean conformation: We obtain the normalized covariance matrix S as for the ENM covariances: Since roto-translational contributions to the variance are approximately removed these matrices effectively have rank 3N − 6, provided sufficient sampling.The relative magnitudes of atomic fluctuations are obtained as where s i 2 denote the atomic fluctuations for atom i and S ii is the 3 × 3 submatrix of S containing the positional covariances of the three Cartesian coordinates for atom i.
B-Factors.When comparing with B-factors we use the isotropic B-factor of each C α atom as they are found in the PDB structures listed in Table 1.These parameters are fitted to experimental observations and typically represent the width of a Gaussian electron density assumed for each atom.The values for this fit are dependent on many factors, like thermal fluctuation, lattice defects, and static disorder.Because of the former they are sometimes taken to be an approximation to the thermal fluctuations of the crystal.We denote the B-factors for a crystal structure by the vector b 2 , with b i 2 the B-factor of atom i. Measures of Similarity.Similarity of Covariances.We benchmark ENMs by comparing the covariance matrices they predict with covariance matrices obtained from MD simulations.We compare such matrices with the rank-normalized Bhattacharyya distance, introduced for this purpose in Fuglebakk et al. 52 For comparing covariance matrices A 1 and A 2 this is where the columns of Q are the q principal components of (A 1 + A 2 )/2 with highest variance and Λ is a diagonal matrix with the corresponding eigenvalues.||denotes the determinant, and we choose q to retain 95% of the total variance of (A 1 + A 2 )/2 with as few principal components as possible.This dimensionality reduction is done to ensure positive definiteness of the matrices.For convenience we report the similarity as the Bhattacharyya coefficient: which is maximally 1, obtained when A 1 = A 2 .Viewing the proteins' Boltzmann distributions as central Gaussian distributions, p A 1 with covariance matrix A 1 and p A 2 with covariance matrix A 2 , eq 24 can be derived from the more general form: See for instance Lee Brestschneider. 53The Gaussian assumption implied in deriving eq 24 from eq 25 is approximate for the MD covariances, while strictly justified for the ENMs.The form used by Merritt 54 to compare anisotropic displacement parameters in crystallographic model refinement is equivalent to eq 24 for Gaussian distributions.For that application, where all matrices are positive definite with the same rank, the rank normalization and dimensionality reduction are not required.
A brief exploration of the sensitivity of the Bhattacharyya coefficient to differences in covariance matrices are provided in the Supporting Information (Figures S1 and S2).
Similarity of Atomic Fluctuations or B-Factors.The atomic fluctuations predicted by MD sampling (eq 22), ENMs (eq 5), or isotropic B-factors can all be expressed as N-dimensional vectors with a value for each C α atom of the protein structure.We compare two such vectors, w 1 and w 2 , by their normalized squared inner product (SIP): Null Model.As a null model for the benchmarking of covariances, the identity matrix would represent the covariance matrix of a protein with uniform fluctuations and no correlations between the atoms.However, the complete absence of correlations can be obtained only by including the rotational and translational degrees of freedom of the whole protein.Therefore we consider instead the projection of the 3N × 3N identity matrix onto the subspace of internal deformations of the protein.The null-matrix, Σ 0 ,is where I is the identity matrix and v are the orthonormal basis vectors for the subspace of global translations and rotations.The matrix is then normalized to its trace as was done for the ENM and MD covariance matrices.For comparing with the GNM we use an analogous definition of the null-matrix, which is the projection of the N × N identity matrix.This is obtained by simply subtracting the outer product of the normalized eigenvector describing uniform displacement, u.The null-matrix for the GNM Z 0 is Note that every entry of uu T is equal to 1/N.These matrices have rank 3N − 6 and N − 1 with all nonzero eigenvalues being the same.Motion represented by these matrices is not completely uncorrelated between atoms, but as the variance is equal along all principal components, collective and local motions have the same energetic cost.The Bhattacharyya measure conveniently allows us to do meaningful comparisons with the null model, even with this degeneracy of eigenvalues. 52D Simulations.For two of the structures we used MD trajectories previously published.We refer to these publications for the modeling details.A 10 ns trajectory of lysozyme (PDB ID: 193L) was obtained with the Amber94 force field 55 by Kneller and Hinsen 56 (the reference reports a 1.2 ns trajectory, which was later extended to 10 ns).A 285 ns trajectory of the GroEL subunit (PDB ID: 1XCK, chain A) was obtained with the Amber03 force field 57 by Skjaerven et al. 58 A 20 ns of trajectory of phospholipase-C (PDB ID: 1PTG) solvated in TIP3P-solvent 59 was obtained with the CHARMM27 force field 60 by co-worker Cedric Grauffel, who kindly offered this trajectory for analysis in this study.Integration was done with NAMD 2.8, 61 with an integration step of 1 fs and a temperature of 300 K. Distance constraints were enforced on all bonds between hydrogens and heavy atoms by the SHAKE algorithm. 62or the remaining four structures we solvated proteins in TIP3P-solvent and simulated using the CHARMM27 force-field.In the ATCase complex zinc ions are coordinated by four proximal cysteine residues, a type of interaction not well represented in CHARMM27.We modeled these interactions with harmonic energy terms restraining the distance between the sulfurs and the ion. 100 ns production runs were obtained for all four systems, using Gromacs version 4.5.4 63with an integration step of 2 fs.Distance constraints on all covalent bonds were enforced by the LINCS algorithm. 64All ensembles were coupled to a thermostat set to 300 K.
Additional details about the MD simulations are provided as Supporting Information.
Analysis of MD Trajectories.For all trajectories we superimposed all frames to the conformation reported in the PDB structure to minimize the RMSD.We then investigated whether this distance was stable across the trajectory.We discarded any leading frames that exhibited very large variation in this parameter.For most proteins no frames needed to be discarded, but for the ATCase-complex and the GroEL subunit as much as 40 and 55 ns were removed, leading to final trajectories of 60 and 230 ns, respectively.Some of the trajectories still had somewhat increasing trends for the RMSD.Nevertheless, we have chosen to use these trajectories, as converged descriptions of collective motion is not attainable.Recent work also indicates that amplitude-normalized ensemble statistics are not very sensitive to simulation time. 32,65To guard from overinterpretation we quantified the heterogeneity of the covariances in a trajectory by comparing covariance matrices obtained from the first and last third of the trajectory.This internal similarity is denoted MD 1/3 .For the analyzed trajectories, plots of the RMSD to the PDB conformation is provided as Supporting Information.

RESULTS
Benchmarking the ENMs.To learn which of the investigated ENMs are better suited to model the inter-residue covariances of compact folded proteins, we calculated the Bhattacharyya coefficient for the comparison of each ENMpredicted covariance matrix (eq 4) to the covariance matrix obtained from MD simulations (eq 21).
The results are tabulated in Table 2. Reported here are also the similarities between the MD covariances and the null model (eq 27) and the similarity between covariance matrices obtained from the first and last third of each trajectory.For the GNM, no result is reported in this table, as the covariance matrix that can be obtained from this model has a different interpretation from those obtained with anisotropic models.For the ANM 0.8 model of the ATCase, the Hessian had rank lower than 3N − 6.These additional zero components were eliminated the same way as the roto-translational components when constructing the covariance matrix (eq 4).
We observe a clear separation between models that have weak long-range interactions (HCA, REACH, ANM 0.8 , and βGM) and models that have relatively strong long-range interactions (pfANM and ANM 1.8 ), the former performing consistently better than the latter.We consider the ANM 1.8 to provide longrange interactions due to the very high connectivity of an atom that interacts with all neighbors within 1.8 nm.Within these two groups there are no consistent trends for ranking the models, although there seem to be a tendency for the refinements introduced in REACH and βGM to have a positive effect.Note that the functional form describing the distance dependence of force constants are varied within these groups, with little effect on the agreement with MD covariance.Similar observations were recently made by Leioatts at al. 36 Since the uniform force constant models are convenient to visualize, we refer to Figure 1 for an illustration of the contrast between weak long-range interaction (ANM 0.8 ) and strong longrange interaction (ANM 1.8 ).
The MD-predicted covariances are dependent on the timescales sampled, the choice of force field, solvent modeling, and sampling parameters as well as important but arbitrary system preparation parameters like initial velocity configurations.Observations regarding the differences between the ENMs that are not consistent for different proteins should therefore be interpreted cautiously.For the case of the GroEL subunit we investigated the dependence of covariances on the initial velocities by comparing 5 trajectories of 100 ns each, all computed from different initial velocity configurations.This was the system with the longest simulation time and the lowest internal consistency (MD 1/3 ).Averaging over all 10 comparisons this gave a mean Bhattacharyya coefficient of 0.84, with a standard deviation of only 0.019.We did not attempt to quantify this dependence for other proteins.Considering the comparison of parts of the trajectories to each other (MD 1/3 ), we note high coefficients indicating that the trajectories give stable estimates of the covariances.
For the GroEL subunit, for which we have data from multiple trajectories, we see that the best performing ENMs are as similar to the MD covariances as covariances obtained from these trajectories are to each other.This is in agreement with earlier observations that a good ENM can be as reliable as a single MD trajectory for capturing the covariance structure of a protein. 16n the other end of the table, we observe that the worst performing models are sometimes only marginally more similar to the MD covariances than the null model.
Alternative Benchmarks: MD fluctuations or B-Factors.We have chosen to benchmark the ENMs by comparing the covariance matrices obtained from coarse-grained ENMs with those obtained from atomistic MD simulations.The primary reason for this is that the covariance structure of the protein determines its collective motions, the kind of motions that ENMs are best motivated for studying.Although the principal components of the covariance matrices are often used to make inference about large amplitude motions, both the ENMs and the MD simulations only explicitly consider small deviations from the equilibrium structure.Under equilibrium dynamics, it is expected that tightly coupled degrees of freedom will remain tightly coupled even on longer time scales, which justifies validating ENMs by assessing the covariance structure.It is however not clear that agreement between atomic fluctuations of models imply agreement between their covariance structures.As other studies have frequently validated ENMs by the atomic fluctuations they predict, we are interested in learning to what extent this benchmark agrees with results for MD covariances.
We therefore benchmarked the ENMs by comparing the atomic fluctuations obtained from the ENMs (eq 5) to both atomic fluctuations obtained from MD trajectories (eq 22) and B-factors.The agreement was quantified using the SIP (eq 26).
Figure 2 summarizes the results and contrasts them with the benchmark against MD covariances.We note that the best models according to the benchmark comparing covariances are also the best in terms of benchmarking with MD atomic fluctuations, except for ANM 0.8 .Despite its high score when benchmarking with covariance, the ANM 0.8 sometimes obtains very low scores both for comparisons with MD atomic fluctuations and B-factors.This is illustrated by the very inconsistent intersect of the green lines with the F and B axes in Figure 2. The ENMs that are in poor agreement with MD covariances are also generally in poor agreement with the atomic fluctuations, while they are in good agreement with the B-factors.Notice in this respect the frequent crossing of red and blue lines connecting C and B axes in Figure 2. Particularly, pfANM has consistently the best agreement with B-factors among the anisotropic models, while almost consistently being among the worst along the other axes.
For benchmarking with MD atomic fluctuations and B-factors, we also obtained results for the isotropic GNM with a cutoff of 0.8 nm.These are shown in gray in Figure 2 (only between the F and B axes), showing that GNM 0.8 quite consistently agrees poorly with the MD atomic fluctuations and agrees well with Bfactors.
Atomic Fluctuations of Short Cutoff ANM.The disagreement between benchmarking with atomic fluctuations and covariances is striking for the case of ANM 0.8 .We inspected the atomic fluctuation profiles for the proteins where this disagreement is larger and found this difference to be explained by the kind of extremely mobile residues that can be expected in under-restrained regions of short cutoff ANM models. 15One example is the GroEL subunit.For this structure ANM 0.8 ranks comparably well when considering covariances but scores poorly both when comparing with B-factors and MD atomic fluctuations (Figure 2).Although most of the profile shapes are in good agreement, there is a very high peak that affects the normalized magnitudes resulting in a low value for their SIP; see Figure S5.This peak consists of three residues in a short loop segment joining two antiparallel β-strands (Figure S6).If the fluctuation of the most mobile residue in this peak is excluded from the atomic fluctuation profiles, the SIP is increased from 0.15 to 0.44.
In the case of the ATCase complex some chain terminals are only connected to two other C α atoms with this cutoff.This leads to unrestrained rotation around a nearby bond, which explains the aforementioned additional zero components of this Hessian; see Figure S7.
Cutoff Dependence for Covariances and B-factors.To investigate the sensitivity of the benchmarks to the choice of cutoff parameter for the ANM, we repeated the benchmark against MD covariances and B-factors for different cutoffs in the range 0.8−2.5 nm. Figure 3 summarizes the results (the large ATCase complex was omitted from this analysis for computational reasons).The cutoff values that are optimal for reproducing B-factors (red marks, Figure 3b) are much higher than those that optimize similarity with the MD covariances (red marks, Figure 3a).Often they are even higher than the 1.8 nm parametrization illustrated in Figure 1.Compared to the MD covariances, trends for B-factors are less consistent between different proteins.Particularly, the differ with respect to how tolerant the reproduction of B-factors will be for lower cutoff values.
The cutoff that optimizes agreement with MD covariances consistently falls in the narrow range of 0.8−1.0nm (Figure 3a).−68 Loss of Collectivity with Long Cutoffs.We observe from Table 2 that the models that agree most with the MD-predicted covariances are those that have weak force constants for the long-  The ATCase complex was excluded from this analysis for computational reasons.
range interactions, and from Figure 2 that stronger force constants improve the agreement with B-factors.Similarly, we see that GNM 0.8 is in agreement with B-factors, but less so with MD atomic fluctuations.
The collective motions of individual protein chains in a crystal are likely to be largely constrained.As acknowledged by the ANM authors, agreement with B-factors (for high cutoffs) is likely to be partly due to these constraints imposed on large-scale structural changes in crystals.Consider, e.g., the difference in constraints on the protein domains between the two different cutoff parametrizations in Figure 1.Similarly we hypothesize that the potential of the GNM (eq 6) constrains the collective motions even for short cutoffs by putting an energetic cost on domain rotation movements.
To investigate the effect of cutoff parametrization on the collective motions, we compared the null model matrices with ANM and GNM covariance matrices, obtained with different cutoffs.Since the null model is representative of nearly uncorrelated motion, with no energy difference between collective and local deformations, we interpret high similarity with the null model as indicative of restricted collective motion.Results are shown in Figure 4.For increasing values of the cutoff parameter the ANM and GNM covariances approach the null model, reaching close to complete agreement at the values obtained from optimizing cutoffs to reproduce B-factors (Figure 3b).
Contrasting the results of GNM and ANM for low cutoffs, we note that the GNM is in close agreement with the null model, even at short cutoffs.For the pfANM, which is the anisotropic model that agrees the most with B-factors and which also have force constants decaying very slowly with distance, comparison with the null model gives Bhattacharyya coefficients in the range [0.97,0.99].Coefficients for the other ENMs are provided as Table S3.

DISCUSSION
The Appropriate Ranges of Interactions.ENMs give an approximate representation of the protein interactions, suitable for investigating thermodynamical properties governed by collective motions in dense protein structures.Validation and benchmarking of such models should therefore make sure to test their ability to represent collective motions.Lacking welldeveloped methods for experimentally determining the covariance structure of proteins, we think that covariances from computational sampling methods are the best standard available for validating ENMs.
Benchmarking the ENMs against covariance matrices obtained by MD simulations, we find that HCA, REACH, ANM 0.8 , and βGM perform consistently better than ANM 1.8 and pfANM.The former four models are in good agreement with MD covariances, while the latter two predict covariance structures close to that of a null model with no collective motion.
The anisotropic ENMs differ only in their force constant parametrization.Since we only treat normalized covariances, the explanation for the differences in benchmarking scores between the ENMs should be sought in the contrast between the force constants within a protein model.The competitive performance of ANM 0.8 reveals that capturing the first coordination shell around internal residues 66−68 already approximates well the covariance structure, even when a uniform force constant is used in this range.This is somewhat surprising given the poor prediction of atomic fluctuation profiles (Figure 2) and the technical problems with low-rank Hessians associated with this parametrization.To more reliably predict atomic fluctuations, some treatment of long-range interactions is necessary.It is important, however, to make sure that there is a proper contrast between the short-and long-range interactions.Simply extending the ranges of the appealingly simple uniform force constant model quickly leads to restrictions on collective motion as illustrated in Figure 4.
Modeling long-range interactions with fast-decaying distancedependent force constants (HCA and REACH) or increasing the granularity of the elastic network (βGM) addresses the problem of hyper-mobile residues seen in the ANM 0.8 .As can be seen from the lack of collective effects with the pfANM, it is necessary to make sure that the long-range interactions decay sufficiently fast.Particularly one should consider if model parameters are robust to variations in protein size and oligomerization, by checking if the system energy is well behaved for lattice sums.Observing this constraint, it seems that the exact form of fastdecaying force constant is not that important, comparing an exponential (REACH) with an inverse power of six (HCA) (see Table 2).Moreover, changing the exponent in a power-decay model like the pfANM increases its performance as the contrast between short-and long-range interactions is increased (see Figure S4).With the exponent equal to −6, performance very close to the HCA is recovered (see Table S4).
The Appropriate Potential.We have indirectly assessed the alternative potential and isotropic assumption invoked by the GNM, by including this model in relevant analysis.It is particularly informative to compare with the ANM, which defines the force constant according to the exact same criteria as the GNM.As the GNM covariance matrix has a different interpretation than those of the other ENMs, we did not make a direct benchmark of covariances akin to the ones found in Table 2.We did, however, compare the atomic fluctuation profiles of the GNM with those obtained from MD and B-factors and found the GNM to behave like the anisotropic models with strong longrange interactions in these respects (ANM 1.8 and pfANM).We also compared the GNM with a null model analogous to what we used for the anisotropic models and found it to have very small collective components.Interpreting this in the light of theoretical considerations, we think that the similarity between the GNM covariances and the null model indicates that the simplicity obtained by the potential in eq 6 comes at the cost of not capturing collective effects to the same extent as the anisotropic models.Figure 4 shows that this is the case even for cutoffs which retain collective motions with the ANM, indicating that it is indeed the potential and not the cutoff that constrains collective motion in this case.
One possible cause for this loss of collective motion is that the GNM puts an energetic cost on rigid-body rotations of domains, in contrast to any physically realistic potential (see for instance Thorpe). 69Many large-amplitude collective motions in proteins contain such rigid-body domain motions.
The Appropriate Benchmarking Standard.In agreement with Leioatts et al. 36 we find that the ENM performance is robust to the choice of formalism for determining force constants.Once a formalism for force constants is chosen, the parametrization depends on the standard chosen for benchmarking.We therefore seek to explain the large differences in collective degrees of freedom between some of the ENMs by contrasting the benchmark against MD simulations with that obtained from Bfactors.
Validating models against B-factors is appealing as these are fitted to experimental observations.The practice can be questioned though, considering neglected effects such as lattice defects and crystal contacts as well as unjustified assumptions made when modeling the B-factors in the structure refinement process.Also, B-factors are are not necessarily proportional to the temperature, as shown by, e.g., Hinsen, 37 indicating that they can not generally be interpreted as thermal fluctuations.More fundamentally, the assumption that the motion of a protein in a crystal is representative of its motion in solution is not justified.Indeed, we see that the ENMs developed to reproduce B-factors have stronger long-range interactions or use a potential that restricts collective motion.This observation has also been reported elsewhere. 34ur analysis shows that very good agreement with B-factors is generally associated with poor agreement with the benchmark against MD covariances (Figure 2).Contrasting benchmarks against MD covariances with benchmarks against B-factors is complicated by the differences between the measures of similarity involved.For a more direct comparison we also compare to the benchmark against MD-predicted atomic fluctuation profiles, using the same similarity measure that was used for the B-factors.These results are in much better agreement with the MD covariances (Figure 2), indicating that the discrepancy between B-factors and MD covariances are due to the quality or applicability of the data source and not the measure of comparison.Moreover, we show that the collective effects in the models that are optimized for reproducing B-factors are so small that the resulting covariance matrix can be very well approximated with the null model (Figure 4).This supports the notion that the strong long-range interactions obtained from optimizing against B-factors reflect the loss of collective motion in the crystalline environment.
An alternative approach to validate ENMs against crystallographic data involves modeling protein crystals as elastic networks. 37The crystal ENM would then not be expected to have larger collective motion than the experimental data it would be benchmarked against.This can possibly be used to parametrize ENMs.The implied assumptions would be that the parametrization is transferrable to a solvated system and that nonthermal contributions to positional variance are nonsystematic.Our current analysis does not address the merits of such methods.
A severe limitation with benchmarking against MD simulations is the intractability of the time-scales at which collective motion can be sampled to its full extent.It has therefore been necessary to analyze trajectories that have not converged in describing the motion that is usually of interest when ENMs are employed.We therefore have to assume that the correlations observed in our MD simulations are similar to those that would be observed in a much longer simulation.Still the discrepancy between the benchmarks obtained from B-factors and MD covariances is consistently related to the collective effects.Since collective effects are a priori expected to be a limitation of the B-factor benchmark and since analysis of the same effects are a main part of the motivation for ENMs, we find the MD covariances to be a more appropriate benchmark.
As the atomic fluctuations of ENMs have frequently been used to characterize the intrinsic dynamics of proteins, it is reassuring to see a very consistent positive correlation between prediction of atomic fluctuations and covariance matrices (Figure 2).For two proteins (phospholipase-C and ATCase), the ENMs that are in poor agreement with covariances give competitive prediction of atomic fluctuations (Figure 2d,g).This serves to illustrate the aforementioned point that agreement with atomic fluctuations does not necessarily imply agreement with covariances.
Sacrificing Simplicity.The simplicity of the ENMs is an appealing aspect, that both aids in interpretation of the results and makes them widely applicable.Their simplicity also implies that there is a large selection of refinements that can be considered to improve the models at the cost of adding complexity.Improvement of coarse-grained models is an active area of research, and the ENMs tested here, as well as others mentioned in the introduction, are just a representative sample of the many variants proposed.As proposed improvements commonly come at the cost of applicability and interpretability, the development of good experimental benchmarks and the careful investigation of generality of models are important.Indeed, from our results we see that even well motivated refinements do not consistently improve the agreement with MD covariances.To be confident that refinements to models actually improve prediction of collective motion, the quality of benchmarks must be carefully scrutinized.

CONCLUSIONS
For analyzing collective motion in proteins, we recommend using models that provide an appropriate contrast between close interactions and long-range interaction.Of the models investigated here, the best examples are HCA, REACH, and βGM.These are able to capture well the covariance structure of the proteins, while still providing reasonable estimates of the atomic fluctuation profiles.
The large difference in performance between the ENMs can be understood in terms of the benchmarking criteria that has been employed in development of the models.In our opinion the discrepancy between the benchmarks obtained against B-factors and MD covariances reflects expected differences in protein mobility between solvated and crystalline conditions.We strongly recommend parametrizing models of isolated proteins

Figure 1 .
Figure 1.Illustration of a uniform force constant ENM (ANM or GNM) of the triple functional domain protein (trio).Left: Cartoon representation of the protein, C α atoms shown as red spheres.Center and right: Elastic interactions (gray bonds) for a cutoff of 0.8 and 1.8 nm, respectively.
3 nm.An approximation to the C β displacements as a function of C α displacements can be derived from this as

Figure 2 .
Figure 2. A comparison between benchmarking with MD covariances (C) and benchmarking with either B-factors (B) or MD atomic fluctuations (F).All three axes extend in the positive direction from the origin.The lines connect scores obtained with the same model.HCA, REACH, and βGM are colored in blue, ANM 0.8 is colored in green, and ANM 1.8 and the pfANM are colored in red, pfANM being the darker red.The gray dotted lines connect scores for the GNM 0.8 , for which no comparison of covariances were done.

Figure 3 .
Figure 3.The dependence of the performance of the ANM for two different benchmarks on the cutoff parameter.(a) The Bhattacharyya coefficients (BC) for comparing the ANM covariances (Σ) and MD covariances (S).The maximal similarity for each protein is shown in red.(b) The SIP between ENM atomic fluctuations (σ 2 ) and B-factors (b 2 ).The ATCase complex was excluded from this analysis for computational reasons.

Figure 4 .
Figure 4.The BC for comparing covariances predicted by the ANM and GNM to the null model, with different values for the cutoff parameter.Black: ANM, BC(Σ,Σ 0 ).Gray: GNM, BC(Z,Z 0 ).

Table 1 .
Proteins Used for Benchmarking

Table 2 .
Bhattacharyya Coefficients Comparing MD Covariances with Covariances Predicted from Different ENMs a The coefficients are maximally 1, the value obtained for identical covariance matrices.The first row, MD 1/3 , shows comparisons between covariances obtained from the first and last third of each trajectory.The last row, NULL, compares MD covariances with the null model.b ANM 0.8 covariance matrix had rank lower than 3N − 6.