Effective Harmonic Potentials: Insights into the Internal Cooperativity and Sequence-Specificity of Protein Dynamics

The proper biological functioning of proteins often relies on the occurrence of coordinated fluctuations around their native structure, or on their ability to perform wider and sometimes highly elaborated motions. Hence, there is considerable interest in the definition of accurate coarse-grained descriptions of protein dynamics, as an alternative to more computationally expensive approaches. In particular, the elastic network model, in which residue motions are subjected to pairwise harmonic potentials, is known to capture essential aspects of conformational dynamics in proteins, but has so far remained mostly phenomenological, and unable to account for the chemical specificities of amino acids. We propose, for the first time, a method to derive residue- and distance-specific effective harmonic potentials from the statistical analysis of an extensive dataset of NMR conformational ensembles. These potentials constitute dynamical counterparts to the mean-force statistical potentials commonly used for static analyses of protein structures. In the context of the elastic network model, they yield a strongly improved description of the cooperative aspects of residue motions, and give the opportunity to systematically explore the influence of sequence details on protein dynamics.


Introduction
Deciphering the motions that underlie many aspects of protein function is a major current challenge in molecular biology, with the potential to generate numerous applications in biomedical research and biotechnology. Although molecular dynamics (MD) hold a prominent position among computational approaches, considerable efforts have been devoted to the development of coarsegrained models of protein dynamics [1]. Besides their ability to follow motions on time scales that are usually not accessible to MD simulations, these models also give the possibility to better understand the general principles that rule the dynamical properties of proteins. 13Å.
The elegant simplicity of the elastic network models (ENM) certainly contributed to their popularity, and they have been successfully exploited in a wide range of ap-plications [2,3]. In these models, the residues are usually represented as single particles and connected to their neighbors by Hookean springs [4,5]. The input structure is assumed to be the equilibrium state, i.e. the global energy minimum of the system. Common variants include the homogeneous ENM, in which springs of equal stiffness connect pairs of residues separated by a distance smaller than a predefined cutoff, and other versions in which the spring stiffness decays as the interresidue distance increases [6][7][8]. In all cases, the equations of motion can be either linearized around equilibrium, to perform a normal mode analysis of the system [9][10][11], or integrated to obtain time-resolved relaxation trajectories [12,13].
Despite their many achievements, purely structural ENM also come with severe limitations. Notably, modeling the possible effects of mutations within this framework usually requires random local perturbations of the spring constants [14], or a more drastic removal of links from the network [15]. A few attempts have been made to include sequence-specificity in the ENM by setting the spring constants proportional to the depth of the energy minima, as estimated by statistical contact potentials [16,17]. However, this approach cannot be extended to distancedependent potentials, for they are not consistent with the ground hypothesis of the ENM, i.e. that all pairwise interaction potentials are at their minimum in the native structure. Other studies have led to the conclusion that the ENM behave as entropic models dominated by structural features, and that the level of coarse-graining is probably too high to incorporate sequence details [5,18]. Still, the chemical nature of residues at key positions can have significant effects on the main dynamical properties of a protein. Hinge motions [19], for instance, obviously require some architectural conditions to be fulfilled, such as the presence of two domains capable of moving relatively independently. But the amplitude and preferred direction of the motion are most likely determined by fine tuning of specific interactions in the hinge region. In proteins subject to domain swapping, the hinge loops have indeed been shown to frequently include residues that are not optimal for stability [20]. The importance of the amino acid sequence has also been repeatedly emphasized by experimental studies of the impact of mutations on the conformational dynamics of proteins [21][22][23].
A major obstacle to the definition of accurate coarsegrained descriptions of protein dynamics lies in the highly cooperative nature of protein motions, which makes it difficult to identify the properties of the individual building blocks independently of the overall architecture of each fold. By condensing the information contained in a multitude of NMR ensembles, we build here a mean protein environment, in which the behavior of residue pairs can be tracked independently of each protein's specific structure. This methodology brings an efficient way of assessing coarse-grained models of protein dynamics and of deriving effective energy functions adapted to these models. In the context of the ENM, we identify a set of spring constants that depend on both the interresidue distances and the chemical nature of amino acids, and that markedly improve the performances of the model.

Results
Dynamical properties of proteins from the perspective of an average pair of residues The mean-square fluctuations of individual residues (MSRF) have been extensively relied on to characterize protein flexibility and to evaluate coarse-grained models of protein dynamics [24], in part because of their widespread availability as crystallographic B-factors. However, since igure 1: Schematic illustration of the apparent stiffness γ. A simple model containing 8 beads connected by elastic springs was subjected to 10 7 integration steps under Gaussian noise. Selected values of < (∆Ri) 2 >, γ and γ • are given in arbitrary units. Individually, the pairs A-B and C-D would be identical, but they experience differently the influence of the other beads. As a result, the C-D pair is effectively more rigid than A-B (γAB < γCD). In both cases, the motions are somewhat correlated, as the apparent stiffness γ is larger than what is expected from the knowledge of their individual motions (γ • ). Beads A and E do not interact directly but the effect of the network on their relative motions is captured by the values of γAE and γ • AE .
the MSRF carry little information about the cooperative nature of residue motions, we propose to examine the dynamical behavior of proteins from the perspective of residue pairs rather than individual residues. Information about the fluctuations of interresidue distances is contained in the data of NMR experiments for numerous proteins, and will be exploited here. We define the apparent stiffness of a pair of residues i,j in a protein p: where k B is the Boltzmann constant, T the temperature, and σ 2 rpij the variance of the distance r between residues i and j, in a structural ensemble representative of the equilibrium state. γ pij is defined up to a multiplicative factor, which corresponds to the temperature. We also introduce the uncorrelated apparent stiffness γ • pij , to quantify the impact of the individual fluctuations of residues i and j on the fluctuations of the distance that separates them. This is achieved by using σ • rpij instead of σ rpij in eq. 1, where σ • rpij is computed after exclusion of all correlations between the motions of residues i and j (Methods).
As illustrated in Figure 1, γ can be quite different from one residue pair to another. Indeed, besides the impact of direct interactions, γ is also strongly dependent on the overall fold of the protein, and on the position of the pair within the structure. To remove the specific influence of each protein's architecture, we define the apparent stiffness in a mean protein environment γ(s, d): Np(s,d) ij (2) where s is one of 210 amino acid pairs, d the discretized equilibrium distance between pairs of residues (d ≤ r pij < d + 0.5Å), M p the number of structures in the equilibrium ensemble of protein p, and N p (s, d) the number of (s, d) residue pairs in protein p. Pairs of consecutive residues were dismissed, so as to consider only non-bonded interactions. The mean protein environment is thus obtained by averaging over a large number of residue pairs in a dataset of P = 1500 different proteins (Methods).
The influence of the distance separating two residues on the cooperativity of their motions can be investigated by considering amino acid types indistinctively in eq. 2. Interestingly, γ(d) follows approximately a power law, with an exponent of about -2.5 (Fig. 2). Finer details include a first maximal value occurring for C α -C α distances between 5 and 5.5Å, i.e. the separation between hydrogen-bonded residues within regular secondary structure elements, and a second around 9Å, which corresponds to indirect, second neighbor, interactions. The high level of cooperativity in residue motions is well illustrated by the comparison of γ(d) and its uncorrelated counterpart γ • (d). Indeed, these two functions would take identical values if the variability of the distance between two residues could be explained solely by the extent of their individual fluctuations. In a mean protein environment, however, γ(d) is about two orders of magnitude larger than γ • (d) at short-range, and the difference remains quite important up to about 30-40 A.
The comparison of γ(d) values extracted from subsets containing exclusively small, large, all-α, or all-β proteins indicates that the content of the dataset has a remarkably limited impact on γ(d) (Supplementary Fig. 1). This distance dependence can thus be seen as a general property of protein structures, a signature of protein cooperativity at the residue pair level. Of course, since γ(d) is representative of a mean protein environment, deviations may occur for individual proteins, according to their specific structural organizations ( Supplementary Fig. 2).
The apparent stiffness γ(s) is computed for each type of amino acid pair s using eq. 2, by considering only residue pairs separated by less than 10Å. As shown in Figure  3A, the chemical nature of the interacting residues is a major determinant of their dynamical behavior. Unsurprisingly, Glycine and Proline appear as the most effective ingredients of flexibility. Pairs involving hydrophobic and aromatic amino acids tend to be considerably more rigid, with γ(s) values up to 6 times larger. These differences originate in part in the individual propensities of different amino acids to be located in more or less flexible regions (e.g. hydrophobic core vs. exposed surface loops). However, there is only a limited agreement between γ(s) and γ • (s) (Fig. 3A-B): the correlation coefficient is equal to 0.71, and γ(s) spans a much wider range of values. Beyond individual amino acid preferences, the specifics of residue-residue interactions play thus a significant role in determining the extent of cooperativity in residue motions.

Accuracy of elastic network models in reproducing the dynamical properties of proteins
The computation of the apparent stiffness of residue pairs in a mean protein environment provides an interesting tool to probe the dynamical properties of proteins. It also generates a very straightforward approach to assess the ability of coarse-grained models to reproduce accurately this general behavior. We focus here on four common variants of the residuebased ENM [25,26], which differ only by the functional form of the spring constants κ. The dependence of κ on the interresidue distance r pij is defined by two parameters: the cutoff distance l d , above which residues i and j are considered disconnected, and the exponent α that determines how fast κ decreases with increasing distances: where H is the Heaviside function. The value of the temperature-related factor a p is obtained, for each protein independently, by fitting the predicted MSRF with the experimental ones. This ensures that the amplitude of the individual fluctuations of the beads in the network is on average equal to that observed in the corresponding NMR ensemble, and that the predicted γ(s, d) values can thus be directly compared with those extracted from the NMR data. We consider the following models: ENM 0 10 , ENM 0 13 , ENM 2 50 , ENM 6 50 . These ENM variants were used to estimate the value of σ 2 rpij for each pair of residues in the 1500 proteins of our NMR dataset (Methods), and to subsequently compute γ(d) and γ(s) from eq. 2.
Strikingly, all ENM variants systematically predict γ(d) values to be lower than the experimental ones, at least up to interresidue distances of 20-30Å (Fig. 2). These models overestimate thus the amplitude of pairwise fluctuations, relatively to the amplitude of individual fluctuations. For example, if two residues in a protein undergo highly correlated motions, the amount of thermal energy necessary to induce a moderate variance on the distance between them will generate high variances on their individual coordinates. Consequently, if the motions of the beads of the ENM are less coordinated, adjusting the scale of the spring constants to reproduce the amplitude of individual fluctuations leads to an overestimated variance on the interresidue distances, and thus to lower γ(d) values. This problem is particularly apparent when κ is assumed to decrease proportionally to the square of the interresidue distance, in the ENM 2 50 . Although this model was shown to perform well in predicting MSRF values [8], our results suggest that it negates almost completely the coordinated aspect of residue motions. Indeed, as shown in Figure 2, the γ(d) values predicted by this model are very close to those obtained from the experimental data after removal of the correlations between the motions of the different residues (γ • (d)). This observation is consistent with the extremely short atom-atom correlation length characteristic of the ENM 2 50 , recently estimated on the basis of an X-ray structure of Staphylococcal nuclease [25].
The ENM is often considered as an entropic model, not detailed enough to include sequence information in a relevant way [5,18]. It is therefore hardly surprising that common ENM variants produce a poor description of the sequence specificities of protein dynamics. Individual amino acid preferences for more or less densely connected regions are responsible for some variety in the predicted values of γ(s) (Fig. 3C). However, this variety is far from matching the one observed in the experimental data, as shown by a much narrower range of γ(s) values, and a limited correlation coefficient with the experimental γ(s) values, e.g. 0.62 for the ENM 6 50 ( Supplementary Fig. 3).

Derivation of effective harmonic potentials
Mean-force statistical potentials are commonly used to perform energetic evaluations of static protein structures [27][28][29]. These potentials do not describe explicitly the "true" physical interactions, but provide effective energies of interaction in a mean protein environment, in the context of a more or less simplified structural representation. Similarly, within the ENM framework, κ(s, d) defines for each pair of residues an harmonic interaction potential. This potential is also effective in nature, accounting implicitly for everything that is not included in the model (e.g. the surrounding water). Hence, we seek to identify the value of κ yielding the most accurate reproduction of the dynamical behavior of each type of pair (s, d) in a mean protein environment, which is conveniently captured by the apparent stiffness γ(s, d).
For that purpose, let us define E bond (s, d) as the energy of the elastic spring connecting two residues of type (s, d), in a mean protein environment: where γ(s, d) is the apparent stiffness extracted from the experimental data. E bond (s, d) is unknown and is expected to be different for different pair types (s, d). The knowledge of γ(s, d) is thus not sufficient to estimate directly κ(s, d). However, from any approximate set of spring constants κ (s, d), we may build the ENM for all proteins in our dataset, to reproduce the mean protein environment, and compute for each pair type an estimated value of the apparent stiffness, γ (s, d), and bond energy, E bond (s, d).
Since the behavior of a given residue pair is highly dependent on its environment, we can make the assumption that E bond (s, d) is a relatively good approximation of E bond (s, d), even if κ (s, d) = κ(s, d): Indeed, if the spring stiffness of a residue pair is underestimated (κ < κ), it will also appear as less rigid in the ENM than in the experimental data (γ < γ). A more detailed discussion is given in Supplementary Note 1. From eqs. 4 and 5, we devise thus an iterative procedure in which κ(s, d) is updated at each step k by confronting the predicted values of the apparent stiffness, γ k (s, d), with the experimental ones, γ(s, d). It is expected to converge when γ k (s, d) → γ(s, d), that is, when the predictions of the model agree with the experimental data: We used this approach to derive, from the NMR data, four novel ENM variants: the distance-dependent dENM ; the sequence-dependent sENM 10 and sENM 13 , with a distance cutoff of 10 and 13Å, respectively, and the sequenceand distance-dependent sdENM (Methods). Interestingly, the κ values for the 210 amino acid pairs in the sENM 10 are relatively well correlated with the corresponding contact potentials [28], even though they result from totally different approaches (Supplementary Fig. 4). Some common general trends can be identified, e.g. hydrophobic contacts tend to be associated with both favorable interaction energies and large κ values (Fig. 4A). However, the overall correspondence remains limited, indicating that the determinants of protein rigidity and stability are related, but distinct. The distance dependence of κ in the dENM is remarkably similar to the r −6 power law that was previously obtained by fitting against a 1.5ns MD trajectory of a C-phycocyanin dimer [6] (Fig. 4B), although our new model presents more detailed features. Notably, κ remains approximately constant up to interresidue distances of 5-6 A, and then drops by about two orders of magnitude to reach a second plateau between 7 and 12Å. The κ values of the sdENM are shown in Figure 4C, for a few amino acid pairs. This model not only combines the strengths of the sENM and the dENM, but also reveals the sequence specificity of the κ distance dependence. The D-R pair, for example, is almost as rigid as I-I at short distances consistent with the formation of a salt bridge, but almost as flexible as G-G at larger distances.

Performances of the new ENM
The sdENM yields a much more accurate reproduction of the dynamical behavior of residue pairs in a mean protein environment than the common ENM variants, as demonstrated by the good agreement between experimental and predicted values of γ(s) (Fig. 5A, Supplementary Fig. 5), and γ(d) (Fig. 5B). Beyond its performances in a mean protein environment, our new model also brings highly notable improvements with respect to previously described ENM variants when it is applied to the specific architecture of a given protein. This is illustrated by two examples, on Figure 5C and Supplementary Figure 6. A more thorough assessment of the ability of the different ENM variants to capture the motions of individual proteins was performed on an independent dataset of 349 proteins. The correlation coefficient between predicted and observed MSRF (r B ) has been widely used in the past but ignores the cooperativity inherent to protein dynamics, and presents other shortcomings. Therefore, we introduce a new measure ( σ ) that quantifies the relative error on the estimation of the variability of the distance between residue pairs, and is thus focused on the cooperative aspects of residue motions (Methods). 50 is better at predicting the individual residue fluctuations (Table 1). Interestingly, the ENM 0 10 , with its simple cutoff distance, appears superior when it comes to the reproduction of cooperative motions ( σ = 0.59). The new ENM variants based on our effective harmonic potentials present enhanced performances in comparison with the common models. In particular, the dENM reaches the same level of quality as the ENM 6 50 for individual fluctuations (r B = 0.69), but surpasses even the ENM 0 10 for the description of cooperativity ( σ = 0.54). On the other hand, the impact of introducing sequence specificity can be examined by comparing sENM 10/13 with ENM 0 10/13 , and sdENM with dENM. It consists in a slight improvement of the correlation coefficient r B , and a pronounced decrease of the error σ , especially at short-(0-15Å) and mid-(15-30Å) range.

Discussion
For the last decades, statistical potentials extracted from datasets of known protein structures [27][28][29] have played a critical role in static analyses of protein structures, with major applications including structure prediction, proteinprotein docking, or rational mutant design. Our study demonstrates that a similar approach can be taken to derive effective energy functions that are specifically adapted to the coarse-grained modeling of protein dynamics.
More precisely, in the context of the ENM, we exploited a dataset of 1500 NMR ensembles to determine the values of the spring constants that describe best the behavior of pairs of residues, as a function of both their chemical na-ture and the distance separating them. The success of our approach is attested by a drastic enhancement of the ability to accurately describe the cooperative nature of residue motions, with respect to previously described ENM variants. Moreover, a definite advantage of our method is that the effective parameters characterizing the strength of the virtual bonds are directly extracted from the experimental data without any a priori conception of their functional form. The fact that the distance dependence of the spring constants that we retrieve is quite similar to the r −6 power law, which was considered so far as underlying one of the best performing ENM variants [6,25], also constitutes a major support to our approach.
In our derivation scheme, the virtual bonds are parametrized so as to reproduce the behavior of amino acid pairs in a mean protein environment. The analysis of the ability of different models of protein dynamics to describe the motions of residues within this environment sheds an interesting new light on the properties of these models. In particular, our results indicate that previous ENM variants underestimate, sometimes dramatically, the rigidity of amino acid pairs at short-and mid-range. Our new model does however provide a much more accurate reproduction of the balance between short-range and longrange coordinated motions. This is arguably a crucial aspect when considering, for example, the consequences of localized alterations induced by ligand binding on signal transduction or global conformational changes, such as in ATP-powered molecular motors.
Importantly, our results also demonstrate that the ENM does not have to be exclusively structural, and that sequence details may be allowed to play a major role in coarse-grained descriptions of protein dynamics. Thereby, this study paves the way towards comparative analyses of motions in proteins that share a similar structure but present differences in sequence. Such investigations will prove particularly interesting in the context of the rational design of (modified) proteins with controlled dynamical properties. Although we focused here on residue-based elastic network models, our approach is not limited to this particular family, and can be readily implemented to evaluate and optimize other coarse-grained models of protein dynamics. Notably, the impact of chemical specificity on the dynamical behavior of residues should be even more accurately rendered by effective potentials based on a more detailed structural description.

NMR Dataset
We retrieved, from the Protein Data Bank [31], ensembles of at least 20 models from solution NMR experiments, corresponding to monomeric proteins of at least 50 residues that present at most 30% sequence identity with one another. Entries under the SCOP classifications "Peptides" or "Membrane and cell surface proteins" were not considered. The presence of ligands, DNA or RNA molecules, chain breaks, non-natural amino acids, and differences in the number of residues per model were also grounds for rejection. These criteria led to the selection of 1849 distinct structural ensembles. A subset of 1500 ensembles was randomly selected for the main analysis, and the remaining 349 were used to assess the performances of the different ENM variants. Unfolded C-or N-terminal tails were automatically identified (MSRF values larger than twice the average for all residues in the protein) and removed from consideration. In each ensemble, the structure with the lowest root mean square deviation from the mean structure, after superposition, is chosen as representative and used to build the ENM.

Elastic network model
The network is built by considering each residue as a single bead, placed at the position of the corresponding C α atom in the input structure, and connecting neighboring beads with Hookean springs [4,5]. The ENM variants considered here differ only by the form of the spring constant κ as a function of interresidue distance and of amino acid types. In all variants, bonded interactions are described by a larger value of κ, defined as ten times the value of κ for non-bonded interactions at a separation of 3.5Å, averaged over all amino acid types. The potential energy of the network is given by: U = i<j (κ ij /2)(r ij − r • ij ) 2 , where r ij and r • ij are the instantaneous and equilibrium distances between residues i and j, respectively. By definition, the input structure corresponds to the global energy minimum, with U = 0. For a protein of n residues, the Hessian H of the system is the 3n × 3n matrix of the second derivatives of U with respect to the spatial coordinates of the residues. The eigenvalue decomposition of H yields the covariance matrix C of the spatial coordinates, which constitutes the output of the model: where the sum is performed over the 3n−6 non-zero eigenvalues λ k of H, and u k are the corresponding eigenvectors. C is a 3n × 3n symmetrical matrix, constituted of n × n submatrices C ij : where ∆x i , ∆y i , and ∆z i correspond to the displacements of residue i from its equilibrium position, along the three Cartesian coordinates. The predicted MSRF of residue i is given by the trace of submatrix C ii :

Variance of the interresidue distance
For each pair of residues in a given protein p, the experimental value of this variance is readily computed from the NMR data: where M p is the number of structures in the NMR ensemble, r pijm the distance between the C α atoms of residues i and j in structure m of protein p, and r pij the average distance over all M p structures. In the context of the ENM, σ 2 rpij values are estimated from the covariance matrix of the spatial coordinates, by standard statistical propagation of uncertainty: where J is the Jacobian of the distance r pij as a function of the six coordinates (x i , y i , z i , x j , y j , z j ). This estimation of σ 2 rpij relies on the validity of the first order Taylor expansion of the distance as function of the coordinates in the vicinity of the average distance. We ensured that no systematic bias arose from this approximation (Supplementary Fig. 7). To quantify the impact of the individual motions of residues on their relative positions, we use eq. 10 to compute (σ • rpij ) 2 in an artificial construct where residue motions are not correlated. This is achieved by extracting the covariance matrix from the NMR data, and setting to zero all submatrices C ij where i = j.

Iterative procedure
The values of the spring constants of the new ENM variants were derived from the dataset of 1500 NMR ensembles using eq 6. For the dENM, sENM 10 and sENM 13 , the initial values of the spring constants were set equal to the experimental values of the apparent stiffness: κ 0 (d) = γ(d) or κ 0 (s) = γ(s). Note that the γ(s) values were computed by considering only residue pairs separated by a distance lower than the cutoff of 10 or 13Å. For the sdENM, the κ 0 (s, d) values were set equal to the final values of the spring constants in the dENM, κ(d), for all amino acid types. A correction for sparse data was devised to ensure that κ(s, d) tends to κ(d) when the number of residue pairs of type (s, d) is too small to obtain relevant estimations of σ 2 r (s, d). Instead of eq. 2, we used the following definition to compute both the experimental and predicted apparent stiffness: N p (s, d) is the number of pairs of type (s, d) in protein p, M p is the number of structures in the NMR ensemble of protein p, and S is an adjustable parameter set to 500.
The κ values were rescaled after each iteration step, so that the average value of κ over all amino acid types is equal to 1 for pairs separated by a distance of 6Å. Residue pairs of a given type (s, d) for which κ(s, d) < 0.001 (after rescaling), were considered to establish no direct interaction: κ(s, d) was set to 0, and they were no longer considered in the iterative procedure. The performances of the new ENM variants after the first nine iteration steps are reported in Supplementary Table 1. The procedure converged rapidly for the dENM and the sdENM, and the final models were selected after 5 and 3 iteration steps, respectively. The sENM variants did not improve significantly with respect to the initial models (k = 0), indicating that κ(s) = γ(s) is a good approximation, contrary to κ(d) = γ(d). The procedure was thus stopped after one iteration step, for both the sENM 10 and the sENM 13 .

Performance measures
The ability of coarse-grained models to accurately describe protein dynamics is commonly evaluated by computing the Pearson correlation coefficient between predicted and experimental MSRF, < (∆R i ) 2 >, over all i = 1, ..., n residues of a given protein: where, for simplicity, B i was used instead of < (∆R i ) 2 >. There is indeed a direct relationship between the MSRF and the cristallographic B-factors: B i = (8π 2 /3) < (∆R i ) 2 >. B exp i and B pre i correspond thus here to the MSRF of residue i extracted from the NMR data and predicted by the ENM, respectively. The scale of the predicted MSRF values depends on the scale of the spring constants, which are only defined up to a constant factor. This factor was determined, for each protein independently, by fitting the scales of the predicted and experimental MSRF, i.e. to ensure that: Although it has been widely used in previous studies, r B is probably not the most adequate measure to evaluate the performances of coarse-grained models of protein dynamics. As pointed out previously [24,25], it does indeed present several shortcomings: e.g. it is strongly affected by the presence of highly flexible regions, and does not account for possible flaws leading to an intercept of the regression line different from zero. Most importantly, the MSRF describe individual fluctuations but provide no information about the cooperative aspects of residue motions. The quality of the MSRF predictions gives thus no guarantee about the ability of the model to describe the cooperativity of protein dynamics. The ENM 2 50 provides an interesting example, for it performs quite well in predicting the MSRF but basically negates all cooperativity (Fig. 2, Table 1).
Therefore, we introduce a new measure that exploits the information contained in the correlation matrix C, to quantify the error on the estimation of the fluctuations of the interresidue distances: (14) where N p is the number of non-bonded residue pairs in protein p, σ (exp) rpij and σ (pre) rpij are the experimental (eq. 9) and predicted (eq. 10) values of σ rpij , respectively. σ (pre) rpij is obtained after fitting the experimental MSRF with the predicted ones (eq. 13). The error is normalized by σ • rpij , which is the expected value of σ rpij given the individual, anisotropic, fluctuations of both residues extracted from the NMR data, but neglecting all correlations between their respective motions. This normalization ensures that the contributions of the different pairs of residues are equivalent, and that the measure is not dominated by highly flexible regions.
Both r B and σ are computed independently for each of the 349 proteins of our test set, and the average values are reported. We also report the short-( SR σ ), mid-( MR σ ), and long-range ( LR σ ) contributions to σ , obtained by considering only pairs separated by 0-15Å, 15-30Å, and more than 30Å, respectively.