A consensus view of protein dynamics

For further details, see SI Data Set and SI Figs. 5–17

Previous to any dynamic study, we need to determine whether force fields are providing a similar picture of protein structure and whether such a picture is similar to that derived experimentally. Analysis of collected samplings indicates an average divergence (α in Eq. 1) of ≈2 Å between force fields (slightly larger deviations are found in G simulations), which, considering the thermal noise of the simulations (Ω in Eq. 2), suggests a similarity of ≈70% between the four samplings and an average “effective distance (Ω⁻¹) between them of only 1.4 Å (slightly larger values are always obtained for G-simulations; see Table 2). This finding indicates that all simulations are, in fact, sampling a similar region of the conformational space. This suggestion is supported by the analysis of the radii of gyration and solvent accessible surface (differences of ≈1% in radii of gyration and 2% in solvent-accessible surface between the four force fields; see Fig. 1).

Not only are the structures sampled in the four simulations similar, but they are also close to the experimental conformations (see Fig. 1). Thus, differences in radii of gyration between MD samplings and experimental structures are <1%, and average differences in solvent accessible surface are ≈5%. The average backbone rmsds between simulated and experimental structures are ≈2.0 Å (A, 1.9; C, 2.0; O, 1.9; G, 2.5 Å), close to the thermal noise of MD simulations. Quite surprisingly, the presence or absence of disulfide bridges does not modify the deviation of our samplings from experimental structures, which is, however, dependent on the origin of the experimental structure. Thus, proteins showing the largest displacements with experiment correspond always to structures solved by NMR spectroscopy (see Fig. 1 and Table 1). The deviations in these cases are due to displacements in regions of large flexibility, where NMR structures are in general poorly defined. This is clear when calculations are repeated using the TM-score, a more robust measure that reduces the weight of flexible entities in the rms-fit. In this case, no system is found for which deviation from experimental data are >2.5 Å (see Fig. 1). Moderate deviations between MD trajectories and experimental model might then reflect either lack of quality in some parts of the experimental model, or just the intrinsic differences between simulation and experimental conditions. Clearly, the commonly made assumption that any deviation from experimental structure of MD simulations is always due to simulation errors is not correct and can ignore important physical characteristics of proteins.

Native protein–protein contacts are well preserved during simulations (in average, 90% correct predictions). The total number of intramolecular hydrogen bonds found in experimental structures is very well preserved (see Fig. 2), but obviously, hydrogen bond donor and acceptors are more promiscuous in MD than inferred from single experimental structures (see below). All MD simulations preserve very well (90%; see Fig. 2) secondary structure (α-helices and β-sheets). In summary: (i) different force fields provide similar picture of protein structure, and (ii) MD simulations sample regions close to experimental structures. Divergences between simulations and experimental structure signal proteins with very flexible moieties whose exact conformation is often not well defined. The two prerequisites noted above are then fulfilled, and we can then safely analyze the dynamics properties of proteins.

B-factors are the simplest method to analyze local deformability and have the advantage to be accessible from x-ray data, which provides us with an additional test on the quality of MD-trajectories. As noted in Fig. 3, there is a very good agreement between experimental and MD-derived B-factors (see Fig. 3). Detailed analysis (see SI Fig. 5) show that such an agreement is especially good for A, C, and O simulations, whereas B-factor distribution is slightly displaced to higher values in G trajectories. In all of the cases, β-sheet segments are the most rigid, whereas turns are in general the most flexible ones, in good agreement to what is deduced from experimental data. In fact, the only significant difference between MD and x-ray B-factors is that, although experimental B-factors above 60 Ų are rare, they are not so uncommon in MD simulations. Inspection of B-factor distribution in the different proteins (see example in Fig. 3 and SI Fig. 10) show a good correlation (average Spearman's rank correlation of 0.7 between experimental and MD distributions). This means that proteins residues that appear as the most flexible ones in MD-simulations are those with the largest mobility in the x-ray structures. It seems that MD-simulation is just allowing the movement of flexible residues that were “frozen” by the crystal lattice. This hypothesis can be tested by computation of crystal effective temperatures (see Eq. 3). Considering harmonic residue oscillations as those with B-factors below 60 Ų (see Fig. 3) we obtain an “effective temperature” of 290 K for the crystal. This demonstrates that the effect of crystal lattice is not, as often assumed, a general “effective cooling” of the general protein. The reduction in kinetic energy (i.e., the reduction of “effective temperature”) induced by the lattice is localized in a few very flexible residues, whose mobility is severely reduced with respect to the situation found in diluted aqueous solution.

Lindemann's disorder index (see Eq. 4) provides a global picture of the flexibility of proteins compared with that of macroscopic solid or liquids (26). Quite independent of the force field and on the presence or absence of disulfide bridges MD simulations suggest that proteins as a whole behave like dense liquids (Δ_L=0.28 ± 0.06). However, this is the combination of the “solid-like” properties of the buried interior Δ_L=0.18 ± 0.04 (0.16 ± 0.04 if only buried backbones are considered) and of pure liquid properties of the exposed side-chains Δ_L=0.38 ± 0.07. In summary, as suggested by experimental studies in particular systems (26, 27), proteins are “melted-solids” with strong differences between the near-solid interior and the full-liquid exterior.

Essential dynamics provides and additional test of the global flexibility of proteins. Interestingly, all of the force fields provide a similar picture of the type of movements that are more important in describing protein similarity (Γ ≈ 85% see Table 2 and SI Figs. 6–15). Space dimensionality (i.e., the number of relevant deformation nodes, ref. 28) depends linearly with the number of residues of the protein (r² > 0.94; see Fig. 4), with no clear changes induced by the presence of disulfide bridges. Given A, C, and O simulations, we can derive a consensus regression equation (dim = 32.5 + 0.58 × Nres; r² = 0.93), which indicates that each residue adds only 0.6 dimensions to the conformational space, indicating that the global structure of the protein strongly limits, in a noncooperative manner, the number of accessible conformations of the constituting peptides.

The first essential deformation movements have associated very low frequencies (<70 cm⁻¹ for the first 25 essential modes in the C_α space) indicating that significant deformations are expected at room temperature along these modes. Entropies derived from these frequencies (see Materials and Methods) depend in all cases linearly with the number of residues (r² > 0.98; see Fig. 4), confirming the lack of cooperative effects in determining the accessible conformational space of mono-domain proteins. Perfect match is found between A, O, and C simulations, whereas G trajectories lead to more disordered structures. Interestingly, the presence of disulfide bridges does not introduce any sizeable effect in determining the entropy of the folded protein. Thus, we can conclude that, in contrast with the general believe, disulfide bridges do not especially restrict the conformational space of native proteins.

A single experimental structure presents a clearly defined pattern of stabilizing interactions (salt bridges, hydrogen bonds, and hydrophobic interactions). Obviously, when flexibility is allowed, such a pattern becomes more diffuse and interacting groups are more promiscuous. Thus, despite the total number of stabilizing interactions preserved in MD, there is, in all cases, a reduction in the number of “permanent” interactions (defined as those occurring for >80% of the trajectory) with respect to those in the experimental structure. Thus, the number of permanent hydrogen bonds and hydrophobic interactions are ≈68 and 81% of those existing in experimental structures. The number of “permanent” salt bridges is ≈77% the experimental ones for A, C, and O simulations, whereas a bigger loss of salt bridges is found in G simulations (see Table 3). Once “permanent” interactions are located, we can determine their contribution to the global stiffness of the protein by deriving the associated force-constant (see Eq. 6). Results in Table 2 show that, in agreement with chemical intuition and irrespective of the force field used, hydrogen bonds are individually stiffer (force constants in the range 22–29 kcal/mol Ų) than hydrophobic contacts (force constants in the range 3–6 kcal/mol Ų). More surprising is the soft nature of salt bridges (29), especially of those involving lysine, which in fact behave similarly (in terms of stiffness) to hydrophobic contacts. MD strongly suggests that individually, hydrogen bonds are the main responsible of protein stiffness, whereas salt bridges (especially with Lys as donor) are quite labile and promiscuous. It is then clear that none of the theoretically more stabilizing interactions disulfide and salt bridges are controlling protein dynamics.

Protein dynamics is complex; even in the hydrophobic core, side chains may undergo rotameric changes on a time scale longer than that studied in the present work, which is limited to study short time scale dynamics. Furthermore, classic force fields are rough approximations to the real potential energy function. However, our massive dynamic analysis provides strong evidence that most of the behavior detected during MD simulations is force field independent and that deviations from experimental structures, when they exist, are not spurious, but signal a physically meaningful process. Proteins behave like biphasic systems, with a solid core surrounded by external liquid layer, whose flexibility is mostly restricted by hydrogen bonds, not by disulfide or salt bridges. At the monodomain level, the complexity of the accessible conformational space increases slowly, but linearly with the size of the protein, with no obvious cooperativity effects. Finally, and quite surprisingly, crystal lattice does not reduce the kinetic energy of all atoms of the protein, but just of a small set of otherwise too-flexible atoms, those that, when the lattice is released, lead to the largest conformational changes.

We built a list of proteins to simulate the 30 most populated folds according to SCOP (30, 31), CATH (32), Dali (33, 34), and Dagget's databases (35) (see SI Data Set). When several good three-dimensional models exist in a class, we favor structures without bound ligands.

As in our MODEL database (http://mmb.pcb.ub.es/MODEL), we used a common, automatic setup procedure designed to guarantee reasonable ionization states, no electrostatic unbalances, and a good hydration before simulation starts. Thus, in all cases, experimental structures (for NMR, the first one deposited in the Protein Data Bank) were titrated to define the major ionic state a neutral pH, neutralized by ions, minimized, heated, and hydrated (with special care in introducing structural waters). These systems were then preequilibrated for 0.5 ns with parm99-AMBER force field, and then equilibrated (0.5 ns) with each force field (see below).

Equilibrated structures were used as starting points for 10-ns production trajectories, performed at constant pressure (1 atm) and temperature (300 K) using standard coupling schemes (the same in all cases). Trajectories for three (Protein Data Bank ID codes 1CQY, 1OPC, and 1KTE) ultrarepresentative structures of α, β, and α/β folds (28) were extended to 100 (A, C, and O) or 50 (G) ns to evaluate the reliability of 10 ns trajectories. Results in SI Fig. 16 strongly suggest that 10 ns is long enough for simulations for many analysis purposes. Particle Mesh Ewald approaches were used to deal with long-range effects (36). Integration of motion equations was performed every 1 fs, the vibrations of bonds involving hydrogen being removed by SHAKE/RATTLE algorithm (37, 38). Four force fields were used in production runs: AMBER (17) (A), CHARMM (18, 19) (C), OPLS/AA (20–23) (O), and GROMOS (24, 25) (G). TIP3P (39, 40) was used as water model for A, C, and O simulations, whereas G calculations were done using the SPC model (41) for coherence with the force field. Simulations were done using parallel versions of AMBER8 (42) (A-simulations), NAMD (43, 44) (C and O simulations), and single processor version of GROMACS (45, 46) (G simulations).

Control simulations were performed to determine whether force fields were able to recognize unfolding conditions and lead to the destabilization of the folded form. For this purpose, 30-ns simulations of 1CQY, 1OPC, and 1KTE were performed using 8 M urea and a temperature of 368 K. OPLS (47) parameters for urea were used for O and A simulations, whereas specific CHARMM parameters (48) were used in C simulations. Due to the lack of specific parameters, these control simulations were not performed using GROMOS force field. Results in SI Fig. 16 confirm that force fields are able to recognize strong denaturing conditions starting the unfolding of the proteins (which might take micro- to millisecond to complete) in the 0- to 30-ns simulation window. The final set with control simulations (10 ns; A, O, C, and G trajectories) was performed for three proteins (Protein Data Bank ID codes 2GB1, 1LYS, and 1UBQ) for which residual dipolar couplings and S² distributions have been reported (see references in the legend of SI Fig. 17) from NMR experiments. Comparison of MD- and NMR-derived measures of flexibility is quite satisfactory (see SI Fig. 17) giving confidence in the quality of our theoretical estimates of flexibility (it is not possible to perform this comparison for representative proteins in SI Table 4 due to the lack of experimental data).

Trajectories were analyzed to obtain structural and dynamic properties. Structural descriptors include backbone rmsd, TM score (49), radii of gyration, solvent accessible surface for all heavy atoms (50), secondary structure (51), intramolecular hydrogen bonds, contact maps, Ramachandran plots (52, 53), solvent contact profiles, and others. The global similarity between the structures collected with the different force fields was obtained by computing the pair-cross rmsd (i.e., the rmsd between all of the snapshots collected in the two trajectories; see Eq. 1), and a related similarity index (see Eq. 2).where N is the number of atoms and M is the number of frames

Local dynamic properties are represented by the B factors, which are compared with x-ray values to obtain a measure of “crystal effective” temperatures (see Eq. 3). Global flexibility was represented by Lindemann's indexes (26) (see Eq. 4) and by essential dynamics protocols considering either Cartesian or mass-weighted covariance matrices (54), which provided us with measures as protein entropy (55–57) or the dimensionality of the conformational space (28). The similarity between the nature of the essential deformation patterns of two trajectories (of the same protein) was determined by Γ index (see Eq. 5; ref. 58), considering a small set of 25 eigenvectors (which represent ≈80–85% of protein variance).where a′ is the most probable nonbonded near-neighbor distance, N is the number of atoms, and Δr_i stands for the fluctuation of atom i from its average position.where n is the number of size of the important space (n = 25 here) and ν stands for one unitary eigenvector. Note that Γ = 1 indicates that both trajectories are sampling the same type of essential movements, whereas Γ = 0 means that they are orthogonal.

The flexibility related to key interactions (hydrogen bonds, hydrophobic contacts, and salt bridges) was computed by assuming the harmonic oscillator model (see Eq. 6). Because the use of harmonic model implies a Gaussian distribution of distances, we computed force constant (K in Eq. 6) for “stable” interactions considered as those found in >80% of the trajectory. A salt bridge was defined when the distance N_ζ(K)/C_ζ (R) − C_γ(D)/C_δ(E) was <6.5 Å, hydrogen bonds were annotated by using standard PTRAJ criteria (42), and hydrophobic contacts were counted when two hydrophobic groups (A, V, P, F, M, I, L, W), which are not neighbors in the sequence, have their C_β atoms separated by <10 Å.where k_b is the Boltzmann constant, T is the absolute temperature, and ΔX is the oscillation in the interaction distance from average values.