it is well established that proteins sample 
 a variety of unfolded 
 and partially folded structures , but there is still substantial uncertainty 
 about the exact nature of the structural fluctuations within folded 
 proteins under native conditions . attempts to probe this conformational 
 heterogeneity have led to a range of conclusions , including ( i ) proteins 
 are surface - molten solids , ( ii ) protein 
 side chains undergo substantial fluctuations in the context of a more 
 rigid backbone , ( iii ) protein backbone motions reorient structural elements and 
 create pockets , and ( iv ) proteins undergo local folding / unfolding transitions . rather than being mutually exclusive , it is likely that all of these 
 forms of heterogeneity are present , albeit on different time scales 
 and to different extents . quantitatively establishing the magnitude 
 and dynamics of these fluctuations is an important goal for advancing 
 our understanding of protein stability and function . here , we 
 combine extensive sampling of detailed molecular models 
 with the theory of markov processes to quantitatively describe proteins 
 conformational heterogeneity on time scales across 8 orders of magnitude . specifically , we use molecular dynamics simulations to generate hundreds 
 of microseconds of time evolution for each of a series of proteins : 
 ubiquitin ( 72 residues ) , rnase h ( 155 residues ) , and -lactamase 
 ( 263 residues , figure 1 ) . results of these 
 simulations serve as input for building markov models , indicating 
 the boundaries between adjacent free energy basins and the rates of 
 hopping between them ( see the methods section 
 for details ) . these markov models provide access to dynamics on tens of microsecond 
 time scales for the largest systems examined and millisecond time 
 scales for the smallest . these time scales are orders of magnitude 
 longer than those accessible by any individual simulation . in addition , 
 our models provide kinetic information that is inaccessible to other 
 enhanced sampling schemes , which typically sacrifice kinetic accuracy 
 by altering simulation parameters  such as the potential energy 
 surface or temperature  to achieve broader conformational sampling . therefore , we can capture processes beyond the reach of previous works and make comparisons with experiments that 
 appropriately account for the time scales a given experimental technique 
 can probe . our atomically detailed simulations also provide information 
 that is not accessible experimentally , such as structural information 
 that is not easily extracted from nmr order parameters and low population 
 states that are beyond the detection limits of crystallography . structures 
 of ubiquitin , rnase h , and -lactamase in ribbon 
 ( top ) and surface ( bottom ) representations . this unprecedented access to long - time relaxation within 
 folded 
 proteins allows us to carefully address questions about the extent 
 and facility of structural rearrangements : how variable are proteins 
 side - chain degrees of freedom ? on what time scales if protein cores 
 were essentially crystalline , we would not expect 
 them to exhibit significant rearrangements on the sub - millisecond 
 time scales that can be accessed with our models and methods . however , 
 an analysis of the distribution of side - chain structures readily reveals 
 substantial variability throughout the proteins examined in this work . for example , the first dihedral angle of every side 
 chain ( called the 1 angle ) in ubiquitin visits at least one 
 alternative rotameric state besides the dominant state seen in the 
 crystallographic structure . we see similar behavior in rnase h and 
 -lactamase ( where 100 and 98% of residues visit alternative 
 1 rotameric states , respectively ) . we find even more variability 
 in the rotameric states of other side - chain dihedral angles , consistent 
 with chemical intuition that the 1 rotamer should be more difficult 
 to sample due to steric constraints from the backbone and the need 
 to displace a larger proportion of the side chain . furthermore , we 
 do not observe unfolding within the time scales captured by our models , 
 so these transitions are occurring in the context of a compact structure . to quantify the degree of structural heterogeneity within each 
 protein specifically , 
 we assign each dihedral angle to the gauche ( + ) , gauche (  ) , 
 or trans rotameric states . we then calculate the shannon entropy ( s ) of a residue aswhere n is the total number 
 of possible conformations ( i.e. , the product of the number of alternative 
 rotameric states for each dihedral angle ) and pi is the probability of conformation i. this entropy will range from zero for side chains that 
 do not visit any alternative rotameric states to 4.4 for a long side 
 chain that spends equal time in every possible combination of rotameric 
 states . calculating entropies for each residue reveals extensive 
 variability 
 throughout each of the proteins examined , including substantial heterogeneity 
 within their cores ( figure 2 ) . we classified 
 residues as part of the core if their solvent accessible surface area 
 is less than 0.1 nm and as part of the surface otherwise 
 ( figure 1 ) . the substantial heterogeneity of 
 core residues is more surprising and leads us to conclude that proteins 
 are more liquid - like than crystalline , even within their cores . as 
 in liquids , this variability is made possible in a dense environment 
 by correlated motions spanning large distances . similar trends have been observed in previous computations 
 and experiments but only for small proteins with low stabilities , 
 such as ubiquitin . the instability of 
 these proteins raises the question of how general such observations 
 are , since more stable proteins could easily have substantially less 
 dynamics . therefore , our ability to examine larger , more stable proteins 
 and observe similar levels of heterogeneity is an important contribution 
 that complements recent work on this subject from a crystallographic 
 perspective . histograms of the entropies 
 of side - chain rotameric states ( left ) 
 and the time scales for transitioning between rotamers ( right ) for 
 surface residues ( top , blue ) and core residues ( bottom , yellow ) . in addition to revealing alternative conformations , our computational 
 models complement experiments by providing insight into the time scales 
 for transitioning between different rotameric states . such kinetic 
 information is currently inaccessible crystallographically and may 
 be challenging to obtain from nmr , since different techniques are 
 often required to capture different time scales , as discussed later . our computational analysis shows that transitions within the proteins 
 cores are often slower than those on the surface ( figure 2 ) . furthermore , some dihedral angles exhibit quite 
 slow transitions ( on the 10 s time scale ) even on the proteins 
 surfaces . by contrast , it is common to assume side - chain dynamics 
 typically occur on ns time scales or not at all . the fact that some 
 of these transition times are approaching the total simulation time 
 also suggests that yet more variability may be present on longer time 
 scales . an important 
 question at this point is whether the apparent disorder suggested 
 by our computer simulations is in fact real . or is the structural 
 heterogeneity in our models the result of errors in the force field 
 used to parametrize the interatomic interactions ? to address 
 this question , we make a quantitative comparison with one of the most 
 direct measures of protein flexibility , namely , nmr order parameters . nmr order parameters are derived from measurements of the autocorrelation 
 function of a unit vector along a particular bond ( relative to the 
 molecule s reference frame ) . specifically , the order parameter 
 ( ) iswhere p2 is the 
 second order legendre polynomial , [t ] is a 
 unit vector along a particular bond at time t , and 
  ...  denotes an equilibrium ensemble average . therefore , one will obtain an nmr order parameter 
 of 1 for a perfectly rigid bond and an order parameter of 0 for a 
 freely rotating bond that loses memory of its past orientation . h bond ) have been measured for 
 many proteins , and order parameters for side - chain methyl groups have 
 been measured for a few proteins . one important consideration 
 in comparing to experiment is the different 
 time scales that different nmr protocols address . for example , the 
 molecular tumbling time sets an upper limit on the time scales accessible 
 to standard relaxation methods . for the systems studied here , therefore , 
 we take care to apply the same time scale limitations when calculating 
 relaxation - based order parameters from our computational models , as 
 described in the methods section . more recent 
 residual dipolar coupling ( rdc ) experiments capture nanosecond to 
 microsecond time scale events , so we use the entirety of our data 
 for calculating rdc - based order parameters from our models . rdc data , 
 however , is currently much less abundant and subject to greater statistical 
 uncertainty ( figure 3 ) . given these limitations , 
 we focus on capturing general trends revealed by rdc measurements , 
 in particular the observation that rdc - based order parameters are 
 generally lower than those determined from relaxation experiments . ( top ) 
 nmr order parameters for ubiquitin s side - chain methyl 
 groups measured from experiment ( filled circles ) and calculated from our models ( open diamonds ) . this perspective 
 highlights that the calculated values are within the range of the 
 experimental values and that the rdc - based measurements are generally 
 lower than the relaxation - based measurements . filled circles are 
 for core residues , and open circles are for surface residues . this 
 perspective better captures the degree of agreement between the two 
 sets of values . we obtain reasonable 
 agreement between calculated and experimental 
 order parameters addressing both of these time scale regimes ( figure 3 ) . for example , the root - mean - square deviation ( rmsd ) 
 between calculated and measured relaxation - based order parameters 
 for ubiquitin is 0.18 and the rmsd between calculated and measured 
 rdc - based order parameters is 0.29 ( pearson s r - values of 0.68 and 0.33 , respectively ) . we judge this level of consistency 
 as reasonable considering ( i ) the simplifying assumptions made in 
 deriving order parameters make exact agreement unlikely , ( ii ) our 
 errors are random rather than being systematically in one direction , 
 ( iii ) the large error in experimental rdc - based measurements ( as described 
 above ) , and ( iv ) we capture the qualitative trend that examining longer 
 time scales reveals more structural heterogeneity . for example , our 
 calculated rdc - based order parameters for ubiquitin are 0.29 less 
 than the relaxation - based order parameters , on average ( compared to 
 a difference of 0.30 in experiment ) . the agreement between computation 
 and experiment is generally greater for core residues , which are also 
 typically less mobile than surface residues . we also predict similar 
 differences between relaxation - based and rdc - based order parameters 
 for the side chains of rnase h and -lactamase ( figure 4 ) , where the average rdc - based order parameter drops 
 by at least 0.10 and 0.06 , respectively , compared to relaxation - based 
 order parameters . nmr order parameters for rnase h and -lactamase s 
 side - chain methyl groups calculated from our models ( open diamonds ) , 
 again showing that the rdc - based measurements are generally lower 
 than the relaxation - based measurements . our results are also 
 consistent with recent crystallographic studies . 
 for example , fraser et al . examined 30 proteins with room - temperature 
 crystallography and found that 37.7% of 
 residues populate an alternative 1 rotamer ( i.e. , not seen 
 in structures solved at cryogenic temperatures ) with a population 
 greater than 20% . we find that 40.1% of residues satisfy this criterion 
 in our sample of three proteins , in reasonable agreement with their 
 findings . importantly , we also capture alternative rotameric states 
 that are less populated and , therefore , are invisible to existing 
 crystallographic methods . just as for side chains , examining longer time scales 
 reveals significantly more structural heterogeneity in the backbone 
 than is detectable on shorter time scales ( figures 5 and 6 ) . for example , the average order 
 parameter for ubiquitin drops from 0.90 for calculations mimicking 
 relaxation - based experiments to 0.70 for calculations mimicking rdc - based 
 experiments . we observe similar drops from 0.89 to at most 0.77 for 
 rnase h and from 0.93 to at most 0.86 for -lactamase . as was 
 observed for side chains , the agreement between computation and experiment 
 for the backbone is generally greater for core residues , which are 
 also typically less mobile than surface residues . nmr order parameters 
 for ubiquitin , rnase h , and -lactamase s 
 backbone n h bonds measured from experiment ( filled circles ) 
 and calculated from our models ( open diamonds ) . relaxation - based order 
 parameters for ubiquitin are from tjandra et al . , and rdc - based parameters are from lakomek et al . relaxation - based order parameters for rnase 
 h and -lactamase are from kroenke et al . and savard et al . nmr order parameters for ubiquitin , rnase h , 
 and -lactamase s 
 backbone n  h bonds from computation and experiment . this is 
 the same data as in figure 5 but better highlights 
 the degree of agreement between the two sets of parameters . filled 
 circles are for core residues , and open circles are for surface residues . the validity of our assertions 
 is supported by reasonable agreement 
 between calculated and experimental order parameters for each system 
 ( figures 5 and 6 ) . for 
 example , the rmsd between calculated and experimental backbone order 
 parameters for ubiquitin is 0.08 for relaxation - based order parameters 
 and 0.20 for rdc - based order parameters ( pearson s r - values of 0.38 and 0.17 , respectively ) . in addition to 
 the factors that limit the agreement between computation and experiment 
 that we described previously , the correlation coefficients for the 
 relaxation - based order parameters are low in part because of the small 
 spread in their values  making the low rmsd a more informative 
 measure of our performance . the experimental order parameters also 
 drop from 0.89 to 0.69 , compared to the computational drop from 0.90 
 to 0.70 , showing that we again capture general trends in the data . we also find reasonable agreement between calculated and experimental 
 relaxation - based order parameters for the backbones of rnase h and 
 -lactamase . the rmsd between calculated and experimental values 
 is 0.08 for rnase h and 0.10 for -lactamase ( pearson s r - values of 0.80 and 0.28 , respectively ) . interestingly , 
 there are two published sets of relaxation - based order parameters 
 for rnase h s backbone from 1995 and 1999 and the agreement between 
 computation and experiment is noticeably better for the newer data 
 set ( rmsd decreases by 0.02 and pearson s r - value increases by 0.08 ) . therefore , it is possible that the agreement 
 between computation and experiment would increase for other proteins 
 if the experimental measurements were repeated with newer , more powerful 
 nmr machines . to our knowledge , rdc - based order parameters have not 
 been measured for rnase h or -lactamase , but we predict that 
 they would be lower than the relaxation - based order parameters by 
 at least 0.1 , as discussed above . we 
 have demonstrated that atomically detailed computational models 
 of proteins employing physically realistic force fields generate dynamics 
 whose scope is more consistent with liquid - like behavior than the 
 crystalline character often attributed to folded proteins . specifically , 
 there is tremendous variability in side - chain conformations , even 
 within protein cores . nmr order parameters calculated from these models 
 are consistent with experimental measurements , supporting the validity 
 of our models . applying these same principles to investigate backbone 
 dynamics leads to the prediction that experiments addressing longer 
 time scales will reveal more extensive heterogeneity in these degrees 
 of freedom as well . we expect methodological advances that further 
 extend the range of accessible time scales to reveal even larger fluctuations , 
 and to enable quantification of their thermodynamic and kinetic consequences . accounting for heterogeneity across these scales will be crucial for 
 deepening our understanding of how proteins function and how we can 
 manipulate protein activity . molecular 
 dynamics simulations were run with 
 gromacs 4.5 on the folding@home distributed 
 computing platform . simulations of ubiquitin 
 were performed as described previously . simulations and markov models of rnase h and -lactamase a total of 1000 simulations 
 were started 
 from pdb 1ubq , for an aggregate of 2.3 ms of dynamics . each simulation was run at 300 k using the amber ff96 force field with the gbsa solvation model . a total of 1000 simulations 
 were started from 1f21 , for an aggregate of 100 s of 
 dynamics . each simulation 
 was run at 300 k using the amber03 force field with explicit tip3p water and 9 chlorine ions to neutralize 
 the charge . a total of 1000 simulations were 
 started from pdb 1jwp , for an aggregate of 81 s of 
 dynamics . each simulation was run at 300 k using the amber03 force 
 field with explicit tip3p water and 7 
 sodium ions to neutralize the charge . all markov models were constructed with msmbuilder . following a standard protocol , every 
 10th conformation from the simulations for each protein were clustered 
 with a k - centers algorithm based on the rmsd between c and c atoms until every cluster had a radius , 
 i.e. , maximum distance between any data point in the cluster and the 
 cluster center  less than 1.2  . then , 10 sweeps of a k - medoids 
 update step were used to center the clusters on the densest regions 
 of conformational space . the remaining 90% of the data was then assigned 
 to these clusters and states with only inbound or outbound transitions 
 discarded . on the basis of their implied time scales , a lag time of 
 10 ns was used for ubiquitin , 20 ns for rnase h , and 2 ns for -lactamase . nmr order parameters ( ) are the plateau value of1where p2 is the 
 second order legendre polynomial ( p2(x ) = ( /2)x  /2 ) , [t ] 
 is a unit vector along a particular bond vector at time t , and  ...  denotes an ensemble average . relaxation - based 
 ( or short time scale ) order parameters were calculated by directly 
 evaluating eq 1 at the molecular tumbling time . rdc - based ( or long time scale ) order parameters were calculated as 
 the long - time limit of eq 1 , as has been done 
 previously . specifically , where x is the x component of the unit vector along 
 the 
 bond of interest and the average is taken across all the conformations 
 sampled ( by using one representative structure from each state in 
 the markov model and weighting its contribution to the ensemble average 
 by its equilibrium population ) . h bond , and 
 side - chain methyl order parameters were measured from a unit vector 
 pointing from the carbon atom to the center of mass of the three hydrogens . all measurements were taken relative to a common molecular reference 
 frame by aligning each conformation to the protein s crystallographic 
 structure . error bars from bootstrapping of our computational models 
 were typically on the order of 0.001 and , therefore , are not shown 
 for visual clarity .
basic 
 principles of statistical mechanics require that proteins 
 sample an ensemble of conformations at any nonzero temperature . 
 however , 
 it is still common to treat the crystallographic structure of a protein 
 as the structure of its native state , largely because 
 high - resolution structural characterization of protein flexibility 
 remains a profound challenge . to assess the typical degree of conformational 
 heterogeneity within folded proteins , we construct markov state models 
 describing the thermodynamics and kinetics of proteins ranging from 
 72 to 263 residues in length . 
 each of these models is built from hundreds 
 of microseconds of atomically detailed molecular dynamics simulations . 
 examination of the side - chain degrees of freedom reveals that almost 
 every residue visits at least two rotameric states over this time 
 frame , with rotamer transition rates spanning a wide range of time 
 scales ( from nanoseconds to tens of microseconds ) . 
 we also report 
 substantial backbone dynamics on time scales longer than are typically 
 addressed by experimental measures of protein flexibility , such as 
 nmr order parameters . 
 finally , we demonstrate that these extensive 
 rearrangements are consistent with nmr and crystallographic data , 
 which supports the validity of our models . 
 altogether , these results 
 depict the interior of proteins not as well - ordered solids , as is 
 often imagined , but instead as dense fluids , which undergo substantial 
 structural fluctuations despite their high packing fraction .