proteins fold into 
 their biologically functional 3d structures 
 by forming cooperative networks of weak interactions that compete 
 against the entropy of the flexible polypeptide chain .
these complex interaction networks hold the key to folding 
 mechanisms and the rational design of 
 new protein folds .
this is because proteins 
 reside in their native or unfolded states for long periods of time 
 ( up to days ) , but the transitions between 
 these two states , and thus the formation and disassembly of the interaction 
 network , seem to occur almost instantaneously . indeed ,
advanced 
 analysis of single - molecule experiments has recently shown that at 
 least some folding transitions take place over periods on the order 
 of a few microseconds , a time scale that is also equivalent 
 to previous empirical estimates of the folding speed limit .
a comprehensive understanding of folding interaction 
 networks would thus entail the characterization of rare events in 
 individual protein molecules at an atomistic level of detail , and 
 with sub - microsecond resolution .
both experiments and simulations 
 have led to significant advances 
 in our understanding of the protein folding process , and their respective capabilities and limitations are such 
 that a combination of the two may lead to insights and cross - validation 
 that could not be obtained using either paradigm alone .
experimental 
 methods can reach atomic - level resolution when investigating the millisecond 
 time scale , but can access sub - microsecond 
 time scales only using coarse - grained spectroscopic probes .
molecular dynamics ( md ) simulations , on the other hand , can generate 
 continuous , atomistically detailed folding and unfolding trajectories , but are computationally demanding , and rely on physical approximations 
 whose range of applicability has not yet been fully ascertained . here
we use a combination of nuclear magnetic resonance ( nmr ) experiments 
 and long - time - scale md simulations to elucidate key elements of the 
 folding process of the single - domain protein gpw .
gpw is a 62-residue 
 entire gene product that folds into an antiparallel + 
 topology in microseconds
. its ultrafast 
 folding  unfolding relaxation rate places gpw in the fast exchange 
 nmr regime over the relevant temperature range , and makes it an attractive 
 target for long md simulations .
a combination of thermodynamic and 
 kinetic criteria suggests this moderate - sized domain folds over a 
 low free energy barrier .
moreover , gpw 
 exhibits distinctly sigmoidal equilibrium thermal unfolding with well - resolved 
 pre- and post - transition baselines that 
 should facilitate the accurate analysis of its unfolding thermodynamics 
 at the atomic level . in experimental 
 studies of the atom - by - atom thermal unfolding behavior of gpw
, we 
 observe a multilayered process in which the large - scale structural 
 changes characterizing the global , two - state - like unfolding transition 
 are superimposed on a more intricate set of atomic- and residue - level 
 structural changes .
we observe a similar level of underlying complexity 
 in coordinated computational studies of the equilibrium unfolding 
 of gpw .
the simulations thus support our experimental results and 
 are also consistent with previous computational studies performed 
 on other proteins .
moreover , from our combined experimental 
 and computational analysis , we infer that the complex structural changes 
 that we observe at the residue level with both methods are intimately 
 connected to the protein interaction network that ultimately determines 
 the folding mechanism .
nmr is a powerful tool for 
 investigating protein conformational changes 
 with atomic resolution .
nmr relaxation dispersion methods , for example , 
 render high - resolution structural information on transient , low - populated 
 folding intermediate states ( i.e. ,  invisible states  ) . in principle
, time resolution limits the application of these methods 
 to proteins with a somewhat slow folding rate ( < 3000 s ) .
recently , this limit has been successfully 
 pushed forward to investigate unfolding fluctuations of native gpw 
 at very low temperature ( 273 k ) , taking advantage of the slow down 
 in gpw folding rate at this temperature .
an alternative approach comes from work on the one - state downhill 
 folding scenario .
one - state downhill 
 folders are single - domain proteins that fold  unfold in microseconds 
 by diffusing down a barrier - less free energy surface at all experimental 
 conditions .
in thermodynamic terms , these 
 domains unfold through a gradual , minimally cooperative unfolding 
 process that results in a broad distribution of 
 structure - specific equilibrium denaturation behaviors .
such remarkable thermodynamic features have 
 been exploited to infer key aspects of the folding interaction network 
 from the cross - correlations among hundreds of atomic unfolding curves 
 obtained by nmr in equilibrium denaturation experiments .
however , one - state downhill folding domains 
 are not widespread , and , in principle , their folding mechanisms could 
 be different from those of other proteins .
the important question 
 is whether the same nmr approach can be extended to the more general 
 case of folding over free energy barriers in which the equilibrium 
 unfolding process is also distinctly more cooperative .
this is an 
 important fundamental question because it has been generally assumed 
 that folding of small proteins follows a two - state mechanism , in which case all probes that effectively report 
 on the integrity of the native structure exhibit the same unfolding 
 behavior and there is no net gain in probing the process with atomic 
 resolution .
in fact , the two - state assumption 
 has been used as justification for interpreting any deviations from 
 the global unfolding behavior that might be observable by nmr as arising 
 from processes unrelated to folding . therefore
, observing unfolding 
 heterogeneity by nmr that is systematically larger than experimental 
 uncertainty , yet fully consistent with the global unfolding process , 
 would provide important experimental evidence that protein folding 
 cooperativity is , in general terms , finite and limited , as predicted 
 by theory and often observed in simulations . in the current study , we performed the 
 nmr analysis of gpw thermal unfolding at ph 3.5 to slow proton amide 
 exchange and thus ensure obtaining full chemical shift ( cs ) assignments 
 when the protein is unfolded at high temperature . at ph 3.5 the thermal 
 denaturation of gpw measured by far - uv circular dichroism ( cd ) , a 
 low - resolution backbone - sensitive technique , is sigmoidal and exhibits 
 a rather sharp transition .
analysis with the simplest two - state model 
 ( ignoring changes in heat capacity , which are small for gpw and do not result in any signs of cold denaturation 
 at this ph ) produces a tm = 329 k and 
 h = 135 kj / mol .
relative 
 to neutral ph , the mildly acidic conditions used here decrease the 
 denaturation midpoint of gpw by about 11 but keep the unfolding 
 enthalpy essentially unchanged ( figure s1 ) .
the slightly lower stability is convenient 
 because it centers the unfolding process within the experimental range 
 ( 273373 k ) and makes the determination of baselines more accurate . to check
whether the mildly acidic conditions 
 induce structural changes in the gpw native state , we determined the 
 3d structure at ph 3.5 ( pdb : 2l6r ) and compared it with that obtained previously at 
 neutral ph .
the comparison between the 
 lowest energy structures for both conditions renders a backbone rmsd 
 of 0.7  ( figure 1a ) , which is comparable 
 to the variability between pairs of structures within the ph 3.5 ensemble 
 ( table s1 ) .
therefore , the structural changes 
 induced by ph are minimal . in addition , the slight destabilization 
 observed by cd seems to arise from an overall increase in net positive 
 charge at this ph , rather than from changes in specific interactions . 

first , overall 
 screening of the electrostatic repulsions at ph 3.5 by addition of 
 1 m salt results in complete recovery of the unfolding behavior 
 observed at neutral ph ( figure s1 ) .
second , 
 the aliphatic carbon cs values for the side - chains of the 11 basic 
 residues in gpw ( 7 r and 4 k ) are essentially identical at both phs . 

the cs insensitivity of the r and k side chains to the protonation 
 status of their acidic counterparts suggests the absence of specific 
 electrostatic interactions in the gpw native structure .
third , we 
 determined the pka values for all the 
 carboxylic residues in gpw and the lone histidine using nmr ( figure s2 ) .
the titration shows that all these 
 residues have pka values ( 4.5  0.14 
 for the five e in the sequence , 3.66  0.15 for the three d , 
 and 6.7 for h ) that are very close to the standard values of the reference 
 amino acids in unstructured peptides .
the lack of significant pka shifts indicates that there is no significant 
 coupling between the charge status of those residues and the folding 
 process of gpw . for the nmr analysis of thermal unfolding we 
 used the n amide , c and c cs for the 
 62 residues in the protein ( a total of 180 atomic unfolding curves ) . 

relative to protons , the cs for these atoms offer far superior performance 
 for the high - resolution analysis of protein unfolding because they 
 are easier to interpret in structural terms , exhibit minimal temperature 
 dependence , and ring - current effects 
 are less significant .
the nmr spectra recorded around the denaturation 
 midpoint ( where changes in populations are maximal ) did not show significant 
 line broadening , consistent with gpw being in the fast nmr exchange 
 regime .
interestingly , the cs denaturation profiles proved to be very 
 heterogeneous ( all cs versus t data are given in table s2 ) .
a large set of the curves ( 102 ) was two - state - like 
 ( 2sl ) , thus showing single transitions , 35 showed two apparent transitions 
 ( three - state - like , 3sl ) , and others showed even more complex patterns 
 ( cp ) . fitting the 2sl curves to a two - state model rendered high variability 
 in tm and h that 
 goes well beyond experimental uncertainty ( (h ) = 62 kj / mol and (tm ) = 9.3 k relative to the 17 kj / mol and 1.6 k that are the respective 
 median fitting errors at 68% confidence ) .
the variability in thermodynamic 
 parameters observed for the 2sl curves results in a distribution of 
 native probabilities that exhibits maximal heterogeneity at the denaturation 
 midpoint .
this observation is qualitatively similar to what was originally 
 reported for the one - state downhill folder bbl . for gpw
the distribution of probabilities is comparatively 
 much sharper , in line with its distinctly sigmoidal global unfolding 
 process .
this result confirms that the atomic heterogeneity is proportional 
 to the overall unfolding broadness of the protein , as has been proposed 
 before , and not an artifact from poorly 
 defined high or low temperature baselines .
three - state fits of the 
 3sl curves also produced high variability ( the parameters for all 
 two - state and three - state fits are shown in tables 
 s3 and s4 ) .
however , the global unfolding obtained as the signal - weighted 
 average of all the cs data is a simple sigmoidal curve that overlaps 
 perfectly with the denaturation profile obtained by a low - resolution 
 backbone - sensitive technique such as cd ( figure 1b ) .
this important control demonstrates that the observed variability 
 in cs represents the true gpw thermal unfolding , which is simple in 
 global terms but inherently complex at the atomistic level .
experimental 
 nmr analysis of the equilibrium thermal unfolding 
 of gpw at atomic resolution .
( a ) superposition of the lowest - energy 
 structure from the nmr ensembles obtained at ph 3.5 ( blue ) and previously 
 determined at ph 6.5 ( cyan ) .
( b ) the global nmr thermal unfolding 
 behavior represented by the second component from the singular value 
 decomposition ( svd ) of the 180 cs ( blue , left scale ) is compared to 
 the unfolding curve measured at low resolution by circular dichroism 
 ( green , right scale ) .
( c ) the 15 different types of atomic unfolding behaviors obtained 
 for gpw from cluster analysis .
clusters 18 have a single transition 
 ( 2sl ) , clusters 913 have two apparent transitions ( 3sl ) , and 
 clusters 14 and 15 include curves with more complex patterns ( cp ) . 

each panel shows a representative experimental cs curve as example 
 ( colored circles ) , the number of cluster elements , and the expected 
 global behavior for reference ( black curve ) .
the latter was calculated 
 fixing the thermodynamic parameters to those of the two - state fit 
 to the blue curve in panel b ( tm = 329 
 k and h = 133 kj / mol ) and fitting the  native  
 and  unfolded  baselines for each probe .
interestingly , the distribution of atomic unfolding 
 behaviors throughout 
 the protein is seemingly random .
for example , grouping the 102 simplest 
 curves ( 2sl ) by atom type only reveals weak trends : backbone - reporting 
 probes (
n and c ) unfold on average 
 more cooperatively ( higher h ) than side - chain - sensitive c , which tend to exhibit less cooperative curves and 
 also slightly higher tm ( table s5 ) . a more detailed analysis using data - clustering 
 tools produced 15 different clusters 
 ( figure 1c ) , each with its own distinctive 
 properties ( tables s6s8 ) .
all these 
 clusters contribute to the overall gpw unfolding signal as they all 
 contain significant numbers of curves ( three or more ) and display 
 comparable changes in cs upon unfolding .
the most characteristic one 
 is cluster 5 , which is largest in size and reproduces the global unfolding 
 behavior closely .
in fact , considered alone , cluster 5 could be construed 
 as solid evidence of two - state equilibrium unfolding behavior .
however , 
 such an argument implies discarding 80% of the curves and 
 71% of the total cs signal .
the elements of the different 
 clusters are spread around the gpw sequence without revealing any 
 obvious structural or topology - based pattern ( figure 2 ) . for example , elements from cluster 5 are found in both 
 -helices , the -hairpin , and the short loops connecting 
 them .
what figure 2 reveals , however , is that 
 cs curves from the same residue are often classified in the same or 
 a similar cluster , indicating that multiple atomic probes provide 
 consistent information on unfolding at the residue level .
another 
 general feature that emerges from figure 2 is 
 the large concentration of cp curves on both tails , which is consistent 
 with looser coupling of the tails to the rest of the protein .
overall , 
 the entire cs data set provides a close look at atomistic complexities 
 of the gpw equilibrium unfolding process that are difficult to obtain 
 using other currently available experimental methods , and which demonstrate 
 that nmr can expose a great deal of structural detail behind the rather 
 sharp equilibrium unfolding process of a barrier - crossing fast - folding 
 protein .
we thus conclude that the atom - by - atom nmr approach is extensible 
 from the one - state downhill folding regime to a more general scenario , 
 confirming experimentally what had been proposed before by theory and providing an excellent opportunity to compare 
 experiments and simulations at the atomic level .
the elements 
 of the 15 clusters obtained from the analysis of the 180 nmr cs curves 
 are overlaid on a schematic representation of the gpw native structure . 

md simulations offer the potential for direct examination of folding 
 and unfolding trajectories at an atomic level of detail . until recently , however , even the longest such 
 simulations did not reach the time scale on which protein folding 
 takes place . in recent years , the development of fast kinetic experiments has led to the discovery of many proteins that 
 approach the s folding speed limit , while special - purpose hardware has extended the reach of continuous , 
 all - atom , explicit - solvent md simulations to the millisecond time 
 scale .
these advances have allowed the 
 execution of single equilibrium simulations encompassing the repeated 
 folding and unfolding of a number of small proteins , shedding light 
 on various aspects of the folding process .
such simulations , however , 
 are based on imperfect physical models ( force fields ) of the interatomic 
 interactions that underlie the dynamics of real biomolecular systems . 
 the ability to compare the results of gpw simulations with corresponding 
 experimental measurements is thus useful not only for developing a 
 deeper understanding of the process of protein folding , but for assessing 
 and improving the accuracy of current force fields .
gpw is indeed a particularly suitable protein model for 
 such purposes , as it has been recently demonstrated by comparing , 
 at the chemical shift level , the properties of a partially unfolded 
 species detected by nmr relaxation dispersion and predicted by coarse - grained simulations . in the work reported here
, we performed a number 
 of md simulations in parallel with the experiments described above . 
 using anton , a special - purpose supercomputer designed for md simulations
, we performed ( i ) a 250 s simulation of 
 the reversible ( un)folding of gpw at 340 k ( the approximate tm in the force field ) to characterize the folding 
 kinetics and ( ii ) four independent 200 s simulated 
 tempering ( st ) simulations to calculate temperature - dependent properties 
 equivalent to those investigated experimentally by nmr . from the comparison 
 between the five runs
, we could assess the simulation convergence 
 and the statistical errors of the calculated quantities .
all simulations 
 were started from an extended conformation , and transitions from disordered 
 conformations to an ensemble of structures consistent with the native 
 state ( c-rmsd below 1  from the nmr structure ) were observed 
 multiple times in each simulation .
visual inspection of the c-rmsd 
 time series for the 340 k simulation ( figure 3a ) reveals a high degree of structural heterogeneity , strongly suggesting 
 that gpw samples several metastable states in addition to the native 
 and unfolded states .
a kinetic - clustering analysis is indeed consistent with the presence of several metastable 
 states that interconvert on the microsecond time scale ( figure 3b ) .
near the denaturation midpoint , the native state 
 appears to be in fast equilibrium with a folding intermediate ( i2 ) 
 in which the two helices and the native hydrophobic core are well 
 formed , but where the -hairpin is largely unstructured .
( a ) 
 time series of the c-rmsd from the first structure of pdb entry 2l6q obtained from a 
 250 s equilibrium md simulation at 340 k started from an extended 
 conformation .
( b ) markov model of the folding free - energy surface 
 obtained from a kinetic - cluster analysis of all the simulation trajectories 
 combined .
the rates of 
 interconversion between clusters , estimated from the equilibrium md 
 simulation at 340 k , are also reported .
( c ) population of the metastable 
 states identified in the cluster analysis as a function of temperature . 

statistical errors for the cluster populations are estimated from 
 comparison of the two independent st simulations .
the cluster populations 
 observed in the equilibrium md simulation at 340 k are reported as 
 circles .
i2 is in slightly slower kinetic 
 exchange with a more unstructured 
 state in which only helix 1 is formed ( i1 ) .
the cluster analysis identifies 
 two additional , partly folded ensembles that are stable on the microsecond 
 time scale : one in which the -hairpin is formed but the helices 
 are disorganized ( m1 ) , and another in which helix 1 is formed and 
 interacts in non - native fashion with a partially formed helix 2 ( m2 ) . 

m1 and m2 can be considered misfolded kinetic traps , and have populations 
 of 15% at 340 k. this model is consistent with a folding mechanism 
 in which helix 1 is formed early on the folding pathway , while the 
 -hairpin is the last structural motif to consolidate . from the st simulations
, we calculated the changes in population 
 of each metastable state as a function of temperature ( figure 3c ) .
the populations at 340 k extracted from st and 
 the constant - temperature md simulation ( circles in figure 3c ) agree to within 10% , with the exception of the 
 unfolded state , which is less populated in the latter .
given the larger 
 amount of data collected , the better convergence properties of st 
 over md simulations , and the good agreement between the different 
 st simulations , this discrepancy ( 1 rt error 
 in relative free energy ) is likely to reflect a larger statistical 
 error in the constant - temperature md simulation . as expected , the 
 population of the fully disordered unfolded state ( u ) increases steadily 
 with temperature , and dominates at the highest temperature . at low 
 temperature , i2 and the native ensemble
are the most populated states , 
 but there is some structural heterogeneity with a substantial population 
 of m2 . i1 and m1 reach maximal population near the denaturation midpoint . 

the relative thermal stability of the different structural clusters 
 reflects their relative order of appearance along the folding pathway , 
 similarly to what has been previously observed for other fast - folding 
 proteins , thus indicating 
 that useful mechanistic information can in this case be obtained from 
 thermal denaturation data .
the structural ensemble obtained from the md simulations 
 at 300 k conforms with the noes measured experimentally at the same 
 temperature ( average violation 0.02  ; only 2 noe restraints 
 are violated by more than 1.5  ) .
the native conformations in 
 the simulations are thus very similar to the experimental nmr structure 
 ( figure 4a ) .
we thus 
 measured the folding - unfolding kinetics of gpw at ph 3.5 using nanosecond 
 infrared t - jump experiments .
the folding relaxation of gpw could be 
 well fit to a single exponential decay for all temperatures . at the 
 same temperature of the standard md simulation ( 340 k )
the relaxation 
 rate measured experimentally is 1/(8.8 s ) ( blue in 
 figure 4b ) , which is essentially identical 
 to the relaxation rate previously measured at neutral ph . for the simulations
, we used the number of helical 
 residues as a proxy for the ir signal and calculated its autocorrelation 
 function over the 250 s equilibrium md simulation ( red in figure 4b ) .
the obtained decay could not be properly fitted 
 to a single exponential decay , but could be well fitted to a double - exponential 
 function with rates of 1/(0.5 s ) and 1/(4.6 s ) .
the 
 fast phase corresponds to helix - melting and -forming events that seem 
 to occur within metastable states in the simulation , and it is thus 
 not associated with folding .
the slow process , which accounts for 
 most of the amplitude and corresponds to the overall folding process 
 in the simulations , is in good agreement with the experimental midpoint 
 relaxation rate measured here at ph 3.5 ( figure 4b ) and in a previous work at ph 6 .
the 
 agreement further improves if we consider the low water viscosity 
 of the tip3p model , which typically results in roughly a 2-fold speedup 
 of the folding kinetics . from a coarse - grained 
 viewpoint
, the md simulations thus appear to provide an excellent 
 description of the global features of gpw folding . comparing long - time - scale 
 simulations of gpw folding with experiments . 

( a ) structural comparison showing the superposition of one representative 
 conformation of the gpw folded state from st simulations ( red ) and 
 the lowest - energy structure from the nmr ensemble at ph 3.5 ( blue ) . 

the experimental ir relaxation kinetics to a 
 final t = 340 k is shown in blue , and its fit to 
 a single exponential function in cyan .
the simulation kinetics is 
 shown in red , its fit to a double exponential function in yellow , 
 and the resulting amplitudes of the fast and slow phases in gray . 

the average atomic unfolding 
 signal is obtained from the singular value decomposition of the 180 
 cs experimental or simulation curves , and then normalized by fitting 
 to a two - state model .
the red circles show the simulation curve and 
 the blue line the nmr experimental curve of figure 1b .
the temperature scale for the simulation ( top ) is compressed 
 by a factor of 2.97 and shifted by 7 degrees relative to the experimental 
 scale ( bottom ) .
the upper bar shows the color - coded specific segments in gpw and the 
 vertical lines the segment tm .
in addition to customary comparisons of native structure 
 and global 
 quantities like folding rate and melting temperature , the combination 
 of the nmr experiments ( figure 1 ) and the st 
 simulations offer the opportunity for a more detailed comparison between 
 simulation and experiment with regard to the thermodynamic and structural 
 properties of gpw as a function of temperature . to achieve this , we 
 calculated the average n nh , c , and c cs for all the conformations observed in the st 
 simulations at each temperature using both camshift and sparta+ .
such average 
 computed chemical shifts are equivalent to the experimental values , 
 which for a protein in the fast nmr exchange regime are also population 
 averages .
the first important observation is that individual computed 
 cs curves are also very heterogeneous . the cluster analysis classified 
 the data in 13 clusters with more than three curves , which display 
 significant differences in cooperativity and denaturation midpoint 
 ( figure s3 ) . in general , the computed data 
 recapitulates the heterogeneity observed experimentally with the exception 
 of the most complex behaviors ( cp curves ) .
figure 
 s3 also shows that the ensemble variability in cs for each 
 temperature is small ( although it increases somewhat at the lowest 
 temperatures due to poorer sampling ) relative to the overall changes 
 resulting from unfolding . as observed before for the nmr experiments , 

however , direct comparison reveals that the 
 simulated unfolding is about 3 times broader than its experimental 
 counterpart , whereas the agreement in denaturation temperature is 
 good ( within 7 k ) .
the much broader melting reflects a substantial 
 underestimation of the unfolding enthalpy in simulations  a 
 phenomenon that appears to be a general ( though
not yet completely 
 characterized ) feature of current biomolecular force fields . when the simulated temperature scale is compressed 
 accordingly , however , the average cs curves overlay almost perfectly 
 ( figure 4c ) . the agreement with experiment 

the computed cs curves can not be fitted to two- 
 or three - state 
 models as accurately as the experiments because they have higher data 
 point fluctuations and a more poorly defined low temperature baseline . 
 as a way to surmount this limitation
, we classified the computed and 
 experimental cs curves by protein segments to compare the average 
 unfolding properties of the various secondary structure elements in 
 gpw . in this analysis
the experimental data shows that the second 
 -strand , which has no tertiary contacts beyond the hairpin , is intrinsically less stable ( midpoint 5 
 k below the global tm ) .
all other segments 
 have tm values closer to the global process , 
 with the first -strand exhibiting the lowest tm among them and the beginning of helix 2 the highest 
 one ( figure 4d ) . at this level ,
most of the 
 atomic heterogeneity is thus averaged out , indicating that the various 
 secondary structure elements unfold almost concertedly , with the exception 
 of strand 2 .
the simulations draw a slightly different picture in 
 which there is more variability among protein segments ( figure 4e ) . the simulations do pick up the low stability 
 of the second strand , but extend it to the whole hairpin region ( 1-2 
 and l2 ) .
moreover , in the simulation helix 1 is more stable than helix 
 2 . to investigate 
 the unfolding of gpw in more depth
this method was originally devised for inferring 
 the interaction network of one - state ( global ) downhill folding proteins 
 from equilibrium nmr unfolding experiments .
the degree of coupling between every pair of residues is thus obtained 
 from the similarities among the unfolding patterns of all the possible 
 cs - curve pairs ( see the methods section for 
 details ) .
a strong coupling between two residues is interpreted as 
 indicative of concerted unfolding , whereas a weaker coupling implies 
 uncorrelated unfolding events ( e.g. , residues with very different tm ) .
primary couplings 
 are between residues interacting in the native structure , and secondary 
 couplings occur when the two residues are indirectly connected through 
 mutual interactions with others .
the approach is directly extensible 
 to the nmr data of gpw because it is atomistically complex ( figure 1 , table s2 ) , and is also 
 transferable to the simulations , since we have computed cs thermal 
 unfolding curves for all the relevant atoms .
figure 5 shows the thermodynamic coupling matrices resulting from 
 the nmr experiments ( left ) and simulations ( right ) .
both the experimental 
 and simulation - based matrices exhibit a complex pattern of high ( red ) , 
 intermediate ( yellow ) , and low ( blue ) couplings , indicating that at 
 the residue level there still is significant variation in unfolding 
 behaviors . as first step in this analysis , we take advantage 
 of the md simulation results to better understand and test the mechanistic 
 connection between the thermodynamic coupling matrix and the underlying 
 folding interaction network .
this is an important issue that could 
 not be examined before because the nmr experimental data does not 
 provide the structural information about the folding mechanism that 
 is readily obtained from the simulated trajectories . along these lines
, 
 we find that the thermodynamic information encoded in the simulated 
 coupling matrix ( figure 5 , right ) agrees with 
 the mechanistic information obtained from the kinetic clustering of 
 the simulations ( figure 3b ) .
kinetic clustering 
 shows that the -hairpin is flexible in the native state , and 
 unstructured in all other clusters ( the only exception is m1 , in which 
 the hairpin is the only formed structure ) .
consistently , the hairpin 
 is virtually uncoupled to everything else ( figure 5 , right ) and has the lowest melting temperature ( figure 4e ) .
thus , the structural element with lowest thermodynamic 
 stability in the simulations forms last during folding and is also 
 the least coupled to other protein regions in the matrix . mapping the 
 folding interaction network of gpw from the residue 
upper panels show the experimental 
 ( left column ) and simulation - based ( right column ) residue 
residue 
 thermodynamic coupling matrix ( highest coupling in red , lowest in 
 dark blue ) , with native contacts identified using hollow black circles . 

the bottom right triangle shows the entire coupling matrix , while 
 the upper left triangle shows only primary couplings ( interacting 
 residues ) .
the structural and sequence distribution of the overall 
 degree of coupling ( defined as the sum of the couplings of a given 
 residue to all other residues ) is shown in the intermediate and bottom 
 panels .
residues were classified in six levels of coupling ( from low 
 to high : dark blue , light blue , cyan , yellow , orange , red ) using data 
 clustering tools .
the two helices can be 
 either formed or frayed in various clusters , 
 but when they are formed , they tend to be formed for their full length 
 ( see i1 , i2 , and m2 in figure 3 ) , in agreement 
 with their strong local coupling ( figure 5 , 
 right ) .
helix 1 , which is better defined structurally than helix 2 
 in all clusters , exhibits stronger secondary couplings .
the presence 
 of interactions between the end of helix 1 and beginning of helix 
 2 in multiple clusters ( native , i2 , and to a lesser extent , i1 ) is 
 consistent with the strong interhelix coupling in the matrix , while 
 the very transient interactions between the protein ends are consistent 
 with their weak coupling observed in the matrix .
residue 
 coupling matrices as an analytical tool for probing protein folding 
 interaction networks .
moreover , the computed coupling matrix provides 
 a picture of the folding process in simulation , in which individual 
 helices can form independently of global folding , that is consistent 
 with the simulated ir relaxation profiles ( figure 4b ) and the kinetic clustering analysis ( figure 3b ) .
this consistency further supports the use of equilibrium 
 denaturation cs measurements as a way to obtain residue 
this exercise is 
 highly instructive from the viewpoint of evaluating and enhancing 
 force field accuracy .
figure 5 shows that the 
 simulated and experimental matrices do indeed share their main features . 
 as in the simulation - based coupling matrix ,
for example , the experimental 
 matrix identifies strong coupling between the end of helix 1 and the 
 loop and the beginning of helix 2 , as well as marginal coupling between 
 the protein ends and between the -hairpin and the rest of the 
 protein .
moreover , high local coupling is found in the helical regions . 
 at the residue level
the similarities are notable to the extent that 
 the simulations correctly capture the polarized distribution of strongly 
 and weakly coupled residues along the length of the two helices and 
 the hairpin ( see the color patterns plotted on the gpw structure ; 
 figure 5 middle ) .
such similarities indicate 
 that the simulations reproduce the major structural and mechanistic 
 aspects of gpw ( un)folding . certain differences between the 
 experimental and simulation - based 
 coupling matrices emerge , however , when inspected in more depth . in 
 the experiment ,
the residues with high overall coupling are more scattered 
 throughout the sequence ( figure 5 , bottom ) 
 whereas in the simulation , coupling tends to concentrate in the two 
 helices ( figure , bottom ) , consistent with 
 the observation of a less cooperative unfolding process ( figure 4d , e ) . in particular
, the simulations do not capture 
 a critical structural node that includes the hinge connecting the 
 hairpin with helix 2 ( l2 , residues 3541 ) and the middle of 
 helix 1 ( h1 m , residues 715 ) .
for instance , l2 and 
 h1 m couple with the other hinge ( residues 1924 ) and the beginning 
 of helix 2 ( residues 4148 ) through primary contacts . strong 
 secondary couplings ( red - orange squares with no contacts ) then link 
 both segments to the remainder of both helices and to some residues 
 in the first half of the hairpin , thus expanding the network over 
 most of the protein . in summary
, while the simulations and the 
 nmr experiments give 
 consistent results with regard to the main aspects of the folding 
 process , with early formation of helices and late formation of the 
 hairpin , the simulated thermal unfolding process is more decoupled 
 at the level of secondary structure elements than that demonstrated 
 experimentally , providing a structural basis for the decreased cooperativity 
 observed in simulations ( figures 3 and 4 ) .
we speculate that the higher stability of isolated 
 secondary structure elements in simulation may result in a free energy 
 surface containing more metastable states than observed experimentally . 

a possible reason for this discrepancy is that current force fields 
 and water models tend to produce unfolded states that are too collapsed , and this collapse may lead to overstabilization of individual secondary 
 structure elements , making the folding transition 
 less cooperative .
theory and simulations 
 have suggested that complex folding mechanisms 
 may still result in global two - state kinetics , but obtaining 
 experimental evidence for this type of conjecture has been challenging , 
 since it requires the use of a method that provides a level of structural 
 detail comparable to that of atomistic simulations .
the nmr experiments 
 reported here constitute a step in that direction , revealing that 
 the apparently simple and cooperative global unfolding behavior of 
 the protein gpw , which is in many ways consistent with two - state folding 
 ( and thus in principle an all - or - none process ) , hides a far richer 
 process at the atomic level .
the combination of experiments and simulations 
 demonstrates that analysis of such atomistic details obtained from 
 nmr data in terms of residue 
residue coupling matrices , a technique 
 that was originally developed for the investigation of one - state downhill 
 folding mechanisms , provides important 
 information related to the folding interaction network and mechanism 
 of gpw .
moreover , the inherent features of gpw folding indicate that 
 the experimental determination of residue  residue coupling 
 matrices by nmr should be broadly applicable to fast folding proteins , 
 even if they fold by crossing a free energy barrier .
practically , 
 the reach of the method will be determined by the interplay between 
 two interrelated parameters .
first , the folding rate of the protein 
 must be considerably faster than the chemical shift time scale ( i.e. , 
 rates higher than 1/(200 s ) ) to ensure fast - exchange 
 nmr conditions .
the second parameter is the intrinsic degree of folding 
 cooperativity , because above a certain threshold the atomic heterogeneity 
 might be indistinguishable from experimental uncertainty . as the rate slows down , folding cooperativity 
 will concomitantly increase , and thus the nmr analysis of sharper , 
 more challenging unfolding transitions will also encounter decreased 
 resolution due to line broadening .
nevertheless , there is nowadays 
 an ample collection of single - domain proteins that have been experimentally 
 identified to fold in microseconds and which contains examples from 
 all major structural classes .
this nmr experimental procedure thus 
 emerges as a new approach to the characterization of folding interaction 
 networks and as a complement to md simulations in addressing the longstanding 
 challenge of characterizing the folding and unfolding of proteins 
 with atomic resolution .
a sequence containing residues 162 
 of the original gpw gene described by davidson and co - workers was subcloned into the expression vector pbat . 

both unlabeled and uniformly n- and c - labeled 
 gpw protein were produced using the same procedures previously described .
thermal unfolding experiments 
 monitored by far - uv circular dichroism were performed on a jasco j-815 
 spectropolarimeter equipped with peltier thermal control using a 1 
 mm path - length cuvette .
the samples were prepared at a protein concentration 
 of 30 m in 20 mm citrate buffer adjusted to ph 3.5 .
nmr samples 
 of c , n - labeled gpw were prepared in a shigemi 
 tube at 1 mm concentration in 20 mm glycine buffer , 0.1 mm nan3 , ph 3.5 , 5% d2o / h2o and 100% d2o . under these conditions
nmr experiments 
 were acquired at 294 k in a bruker avance iii 600 mhz spectrometer 
 equipped with a triple resonance z - axial - gradient probe .
sequence 
 backbone chemical shift assignments were obtained from the following 
 experiments : [ h  n ] hsqc , hnco , hncacb , 
 and cbca(co)nh .
side chain h , n , and c assignments were obtained from hbha(co)nh , h(cco)nh , c(co)nh , 
 and hcch - tocsy .
noe data were obtained from 3d n - hsqc - noesy 
 and 4d c - hmqc - noesy using a mixing time of 110 ms .
all 
 experiments were processed with nmrpipe and analyzed with pipp . for the structure 
 determination of gpw , experimental restraints derived from noe cross - peaks , 
 hydrogen bond distances , and dihedral angles were obtained . only unambiguous 
 interproton distance restraints
were used , and errors of 25% of the 
 distances were applied to obtain lower and upper distance limits . 

hydrogen bond distance restraints ( rnh  o = 1.92.5  , rn  o = 2.83.4
 ) were defined according to the experimentally 
 determined secondary structure of the protein .
the program talos+ 
 was used to obtain 94  and  backbone torsion angle constraints 
 for those residues with statistically significant predictions .
structure calculations were performed with the 
 program x - plor - nih 2.16.0 by minimizing a target function that includes 
 an harmonic potential for experimental distance restraints , a quadratic 
 van der waals repulsion term for the nonbonded contacts , and a square 
 potential for torsion angles .
starting 
 structures were calculated and heated to 3000 k and cooled in 30  000 
 steps of 0.002 ps during simulated annealing .
the final ensemble of 
 20 nmr structures was selected based on lowest energy and no restraint - violation 
 criteria .
these conformers had no distance restraint violations and 
 no dihedral angle violations greater than 0.3  and 5 , 
 respectively .
structures were validated using procheck - nmr and molprobity , which show that the family of 20 structures are of considerable 
 high quality in terms of geometry ( 92.6% of the residues populating 
 the most favored regions of the ramachandran plot for the whole ensemble ) 
 and side chain packing .
coordinates were deposited in the protein 
 data bank with accession code 2l6r and chemical shifts were deposited in 
 the biological magnetic resonance bank ( bmrb code 17322 ) .
we obtained the chemical shift values 
 of the side chain co carbons of all asp and glu residues in gpw by 
 analyzing a series of h(ca)co spectra , a 2d version of the hcaco experiment .
the h1 chemical shifts of 
 the his15 were determined from a series of [ h  c ] hsqc spectra .
all the spectra 
 were recorded on an 1 mm c , n - double 
 labeled gpw sample , 0.1 mm dss in d2o at different phs 
 ranging from 1.5 to 8.1 ( 8.11 , 7.16 , 6.01 , 5.09 , 4.08 , 3.57 , 3.02 , 
 1.98 and 1.53 ) and 294 k. nmr spectra were analyzed with sparky .
the pka values 
 were obtained by a nonlinear least - squares fit of the experimental 
 ph titration curves to the following equation : where 1 and 2 are 
 the chemical shift values at the lowest and highest ph , 
 respectively .
data fitting was done with origin pro 8.0 ( originlab ) . as a first 
 step ,
the reversibility of gpw thermal unfolding at the protein concentrations 
 required for multidimensional nmr experiments was assessed using a 
 nmr sample consisting of 1 mm n - labeled gpw in 20 mm 
 glycine buffer , 0.1 mm nan3 , ph 3.5 , 5% d2o / h2o .
this sample was heated up to 100 c for a period of 
 2 h and then cooled back down . the perfect cross - peaks superimposition 
 of the two [ h  n ] hsqc spectra recorded 
 before and after heating , indicated full reversibility in the relevant 
 conditions .
equilibrium thermal unfolding was investigated using multidimensional 
 nmr experiments to monitor the changes in h , n , and c chemical shifts in the 273371 k temperature 
 range .
for this purpose we used a bruker 5 mm txi - probe with z - axis gradient able to stand up to 423 k and a bto temperature 
 control unit .
accurate temperature calibration was carried out monitoring 
 the ethylene glycol cs at every experimental temperature using a sample 
 of ethylene glycol prepared in a regular 5 mm diameter nmr tube and 
 using an air flow - rate of 535 l / h and following the procedure described 
 by amman et al .
briefly , the neat ethylene 
 glycol sample was equilibrated for 30 min at each temperature before 
 acquiring the spectrum .
standard one - pulse 1d h nmr experiments 
 were used to monitor the changes in proton frequency of ch2 and oh with temperature .
these changes were 
 converted to temperature using equation : t(k ) = 466.5 
 
102.00 . for the thermal unfolding experiments 
 we employed a nmr sample at 1 mm concentration of c , n - labeled gpw in 20 mm glycine buffer , 0.1 mm nan3 , 0.01 dss , ph 3.5 , 5% d2o / h2o .
in these experiments 
 the water resonance was used as spectral reference since the observation 
 of the standard dss signal is precluded in multidimensional nmr experiments 
 because of isotope filtering .
thus , before recording the set of multidimensional 
 experiments , a 1d h nmr experiment was recorded and employed 
 to calibrate the position of the water signal relative to dss at each 
 temperature .
pressurizable tubes manufactured by new era enterprises 
 ( u.s . ) with 0.77 mm thick glass were used to minimize evaporation 
 of the sample .
we measured the ph of nmr samples at the different 
 experimental temperatures using a thermal block and a ph - meter recalibrated 
 at each temperature .
we observed that the ph of the samples was held 
 constant at a value of 3.5  0.1 over the whole temperature range . 

protein backbone amide n , h , c , and c chemical shifts were 
 assigned with the following set of experiments : [ h  n ] hsqc , cbca(co)nh , and hncacb .
gpw folds in fast exchange 
 regime relative to the nmr chemical shift time scale , and thus , each 
 atom produces a single chemical shift corresponding to the dynamically 
 averaged conformational ensemble at every condition .
the ensemble 
 average chemical shifts for all relevant atoms were fully assigned 
 in the 273371 k temperature range at fixed intervals of 5 
 k. at temperature conditions within the unfolding transition region 
 ( slightly below and above the tm ) the 
 temperature interval for nmr experiments was reduced to 3 k to increase 
 the data density in the region of largest change in cs .
all experiments 
 were processed with nmrpipe and analyzed with pipp to obtain accurate 
 chemical shift determinations using a multilevel contour averaging 
 procedure . at ph 3.5 and high temperatures gpw experienced slow proteolysis 
 at specific positions in the primary sequence involving glutamine 
 residues .
we thus monitored protein integrity by comparing [ h  n ] hsqc spectra recorded before and after every 
 set of triple - resonance nmr experiments . at temperatures higher than 
 345
k the gpw proteolysis rate was higher than 1/(5 days ) so we used 
 a freshly prepared sample for recording the suite of triple resonance 
 nmr experiments at each temperature above 345 k. folding - unfolding 
 relaxation kinetics of gpw was measured using a 
 nanosecond laser - induced t - jump apparatus equipped with infrared detection . 

samples were prepared at 4 mg / ml concentration in 20 mm deuterated 
 glycine buffer , after multiple cycles of liophilization and dilution 
 in d2o to achieve complete deuteration of the exchangeable amide protons . 

the buffer and sample were adjusted to pd = 3.5 , correcting for the 
 isotope effect on the glass electrode readout .
the relaxation kinetics 
 of gpw after a jump of 10 k to a final temperature of 328 
 k were recorded on a custom - built infrared laser - induced temperature 
 apparatus described previously .
briefly , 
 the fundamental wavelength of a nd : yag laser ( continuum surelite i ) 
 operating at a repetition rate of 4 hz is shifted , by passing through 
 a 1 m path raman cell ( lightage inc . ) filled with a high pressure mixture 
 of h2 and ar , to 1907 nm , resonant with the infrared absorption 
 of d2o . in this way pulses
< 10 ns with energy of up 
 to 30 mj were obtained , that induced a local temperature jump of about 
 1012 when focused onto the sample .
the folding - unfolding 
 relaxation kinetics of gpw was observed using a quantum cascade laser 
 ( daylight solutions ) tuned to 1632 cm to match 
 the absorption of the  helix signal within the gpw amide i 
 band .
the light transmitted from the sample was recorded using a fast 
 mct detector ( kolmar technologies ) coupled with an oscilloscope ( tektronix 
 dpo4032 ) .
samples were held in a cell formed by two mgf2 windows separated by a 50 m teflon spacer and thermostated 
 at the proper base temperature using two peltier thermoelectric coolers 
 ( te technology inc . ) in a custom - built sample holder .
the transmission 
 of d2o was used as an internal thermometer to measure the 
 amplitude of the temperature jump .
molecular dynamics ( md ) simulations 
 of gpw were initiated in an extended conformation .
the 62-amino acids 
 sequence mvrqeelaaaraalhdlmtgkrvatvqkdgrrveftatsvsdlkkyiaelevqtgmtqrrrg 
 was solvated in a cubic 60  60  60  box , 
 containing 6345 water molecules and 15 mm nacl .
asp , glu , arg , lys , 
 and his residues were treated as charged .
the charmm - h force field ( see som ) was used to represent the system .
the system was initially 
 equilibrated at 340 k and 1 bar in the npt ensemble .
production runs 
 were performed with the anton specialized hardware using a 2.5 fs time step .
bonds involving hydrogen atoms 
 were restrained to their equilibrium lengths using the m - shake algorithm .
nonbonded interactions were truncated at 10 
  , and an ewald method was used to describe the long - range electrostatic 
 interactions .
four simulated tempering simulations of 200 s each , were performed 
 in the npt ensemble .
two simulations used 18 temperature intervals between 
 278 and 398 k , and two used 21 temperature intervals between 278 and 
 438 k. exchanges between temperatures were attempted every 20 ps of 
 simulation . to gain additional kinetic information
, we also performed 
 a 258 s equilibrium md simulation at 340 k in the nvt ensemble . 

a kinetic clustering analysis was used to identify the most relevant 
 metastable states on the free - energy surface through a fit of the 
 autocorrelation functions of 400 cc contacts 
 on time scales ranging from 0.5 to 50 s . in order to obtain consistent clustering ,
helical residues for 
 the analysis of the autocorrelation function were identified according 
 to the stride definition , as implemented 
 in vmd . as a first step , 

all obtained experimental cs curves were classified in three groups 
 according to the unfolding behavior : ( 1 ) 102 curves were classified 
 as 2-state - like ( 2sl ) curves because they had a distinctly sigmoidal 
 trend characterized by a single transition and a single peak in their 
 derivative ; ( 2 ) 35 curves were classified as 3-state - like ( 3sl ) because 
 they exhibited two apparent transitions and thus a bisigmoidal shape 
 and two peaks in their derivative ; and ( 3 ) 43 curves did not fit any 
 of the other two categories ( whether because they exhibited multiple 
 transitions or no obvious transition ) , and were thus included in the 
 complex - pattern ( cp ) group .
the cs curves classified as 2sl and 3sl 
 were fitted to a thermodynamic two - state or a three - state unfolding 
 model equation , respectively , using the procedures described previously .
the 
 data set containing 
 all of the experimental cs curves was analyzed globally using data 
 clustering methods .
cs curves belonging to the 2sl and 3sl groups 
 were clustered according to both the thermodynamic parameters obtained 
 in their fits ( tm and h for 2sl and tm1 , tm2 , h1 , and h2 for 3sl ) or according to the probability of 
 the native state as a function of temperature obtained from the two 
 and three fits . both methods rendered equivalent results .
cs curves 
 included in the cp group were clusterized according to their similarity 
 after normalization of the signal using the z - score procedure .
all clustering routines were performed with 
 the kmeans algorithm by running 10  000 trials to guarantee 
 convergence and selecting the solution with minimum residuals .
the 
 102 2sl curves were clustered into 12 clusters , out of which 8 clusters 
 had at least 3 elements , whereas the other 4 clusters had only 1 or 
 2 elements .
clustering of the 35 3sl curves produced 5 clusters with 
 at least 3 elements .
cp curves were grouped in 10 clusters by kmeans 
 according to their zscored shape , out of which only 2 clusters had 
 3 elements or more . using these procedures , a total of 15 characteristic 
 clusters were identified with at least three elements containing in 
 total 132 atomic unfolding curves .
the thermodynamic coupling index 
 ( tci ) for each pair of residues in a protein is calculated from the 
 pairwise comparison between all the cs curves of the first residue 
 in the pair with all the curves of the second residue .
the thermodynamic 
 coupling index of residues x and y is then calculated by summing all the possible pairwise comparisons 
 of atomic folding curves as reported in the following equation : where pm and pn are row 
 vectors from a matrix including all the zscored cs curves in the protein . 

the numerator thus corresponds to the mean rmsd for all cs curves 
 in the data set .
px , i and py , j are the vectors 
 from the same matrix that correspond to the cs curves of residues x and y , respectively ( i runs over all curves of residue x and j over all curves of residue y ) .
the tci for residues x and y is positive when the average euclidean 
 distance between all their cross - pairs is smaller than the mean rmsd 
 for all cs curves in the protein and negative otherwise .
the tci matrix 
 is constructed by repeating the same procedure over all possible residue 
 pairs .
this procedure was applied to both the experimental cs curves 
 measured by nmr and the synthetic cs curves computed with either camshift or sparta+ from 
 all the conformations observed at each temperature in the st simulations .
 the 
 integration of atomic - resolution experimental and computational 
 methods offers the potential for elucidating key aspects of protein 
 folding that are not revealed by either approach alone . here , we combine 
 equilibrium nmr measurements of thermal unfolding and long molecular 
 dynamics simulations to investigate the folding of gpw , a protein 
 with two - state - like , fast folding dynamics and cooperative equilibrium 
 unfolding behavior . 
 experiments and simulations expose a remarkably 
 complex pattern of structural changes that occur at the atomic level 
 and from which the detailed network of residue  residue couplings 
 associated with cooperative folding emerges . such thermodynamic residue  
 residue 
 couplings appear to be linked to the order of mechanistically significant 
 events that take place during the folding process . 
 our results on 
 gpw indicate that the methods employed in this study are likely to 
 prove broadly applicable to the fine analysis of folding mechanisms 
 in fast folding proteins . 
Introduction
Results and Discussion
Conclusions
Methods