the energy landscape theory of protein folding(1 ) asserts that proteins evolved to minimize frustration during the folding process.(2 ) this suggested that the energy landscape of proteins is globally funneled toward the native state , so that many routes are available for the folding polypeptide chain to reach the native conformation . thus , the folding process involves the formation of ensembles of increasingly ordered structures as it proceeds down the funnel . characterization of the main conformers of the protein and the ways taken to move from one to the other is a fundamental objective of protein folding research . while experiments provide the primary evidence regarding protein folding and its determinants , molecular dynamics ( md ) simulations can potentially add up important information thanks to the atomistic resolution they provide.(5 ) ideally , if the force field used were accurate enough and if computational time was sufficient , md simulations could directly identify all the protein conformers and transitions between them . even in that case , however , we would be challenged by the enormous amount of data which constitute the md trajectory . how can we extract the relevant information from the trajectory and provide it in a a tractable representation of the dynamics of the protein corresponds to providing , e.g. , populations of conformers and their rates of exchange . a large variety of methods on one hand , there are approaches based on the projection of the dynamics along a series of observables ( i.e. , the end - to - end distance , the radius of gyration , root mean squared deviation ( rmsd ) from reference structure , dihedral angles , etc . ) . such methods have the advantage of providing the results with a clear and immediate meaning , but may lead to missing the relevant degrees of freedom , as we project an intrinsically multidimensional system onto the few dimensions represented by the observables . on the other hand , another class of approaches is based on a discretization of the trajectory via clustering , where conformations of the peptide chain are grouped together when they are similar according to some metrics defined directly in the multidimensional space of the peptide chain . these approaches reduce the risk of missing relevant degrees of freedom and allow for the determination of purely kinetic properties such as commitment probabilities and mean first passage times . however , they usually require further complex steps of analysis like markov state models or procedures to recover the free energy basins to reach an accurate and concise description of the system . markov state models ( msm ) have received increasing attention as a powerful tool to analyze and integrate md trajectories and obtain a concise and unbiased view of the dynamics of biopolymers . on the basis of the assumption that on a long enough time scale the biopolymer dynamics can be represented as a markov process , these procedures naturally provide a way to identify the relevant macrostates of the system and their exchange rates . the procedures usually require an initial discretization of the trajectory into microstates that have to satisfy certain criteria ; e.g. , conformations separated by free energy barriers have to belong to different microstates,(13 ) and the size of the microstates has to be sufficient to provide a statistically significant estimate of the transition probabilities into the other microstates . although msms provide tools to check a posteriori the markovianity of the system , other ways to cross - check the results and identify possible weaknesses in the model can improve the overall description and understanding of the biopolymer dynamics . here , we analyze extensive md simulations of a model -hairpin , for which hundreds of folding and unfolding events could be sampled , thanks to the adoption of a continuous model for the solvent . using a bayesian approach , we optimize a msm to fit the trajectories projected onto two relevant observables , and from the eigenvectors of the rate matrices we extract information on the macrostates involved in the slowest relaxation dynamics . , we carry out a kinetic network analysis which consists of building a network ( or graph ) out of the clustered trajectories by associating the conformational clusters with the nodes of the network and the observed transitions with the links between nodes . we then merge the results of the markov state analysis and those of the kinetic network analysis . this is done by projecting the markov macrostates on the network representation and using the markov macrostates as target states for folding commitment probability calculations on the network . we show that this combined approach allows for the determination of the main conformers of the peptide , their exchange rates , and also the transition state ensemble , and helps to provide a concise view of the peptide dynamics . the methods described are found to be particularly useful in the present case , as they allow for the identification of different conformers within the denatured state ensemble . this is especially important because partially unfolded and denatured states of proteins have been implicated in several cellular processes , and the importance in cellular function of natively unfolded and marginally stable proteins is being increasingly recognized.(23 ) similarly , the role played by the residual structure of the denatured states in determining the folding process needs to be carefully evaluated.(24 ) we finally provide a structural and thermodynamic description of the states identified for the hairpin peptide . the 12-residue -hairpin tryptophan zipper 1 ( tz1 , figure 1 , pdb i d 1le0)(25 ) with a type ii turn and a core of four interlocking tryptophan residues was modeled using implicit - solvent molecular dynamics simulations . simulations were run using the program charmm,(26 ) with all heavy atoms and polar hydrogen atoms being modeled explicitly using a united - atom parametrization ( param19 ) . the aqueous effects of the solvent on the solute are described by a mean - field approximation based on the solvent - accessible surface area ( sasa ) of the peptide.(27 ) to improve sampling , solute  solvent frictional forces have been neglected ; this has no appreciable effect on the thermodynamic properties of the system . the 12-residue -hairpin tz1 : from the protein data bank with ( a ) backbone carbon ( cyan ) , nitrogen ( blue ) , oxygen ( red ) , and hydrogen ( white ) atoms and residue labels , and ( b ) tryptophan side - chains and secondary structure , with gray for random coil , red for extended -sheet , and tan for turn sequence . ten md simulation trajectories of 2 s each and different initial velocities were run at 300 , 330 , and 360 k. the temperature was kept constant during the simulations by weak coupling to an external bath , as described by berendsen et al.,(28 ) with no pressure constraints being imposed on the system . the shake algorithm(29 ) was used to fix covalent bond lengths , allowing a 2 fs integration time step . the trajectories were examined by observing the deviation from the nmr structure in the root - mean - square distance ( rmsd ) between backbone c atoms . the number of backbone hydrogen bonds ( h - bond ) was also examined for the peptide . instead of using a step function to describe a h - bond as either formed ( 1 ) with an interatomic distance of less than 2.6  or not formed ( 0 ) with an interatomic distance in excess of this cutoff , a sigmoidal function of the distance d was used fhb(d ) = 1/(1 + ( x/2.6 ) ) . using this method , a continuous range of values between 1 and 0 is obtained , with fhb(d = 2.6  ) = 0.5 representing a conformation with an interatomic distance of exactly 2.6  , fhb(d ) values close to 1 representing conformations with interatomic distances less than 2.6  , and values close to 0 representing conformations with interatomic distances greater than 2.6  . five native h - bond contacts were identified along the backbone of the native conformation . since the h - bond nearest to the turn sequence was formed for conformations in all states , the remaining four h - bonds were taken as a reference by which to assess the degree of nativeness of a given structure , through the variable hb = ifhb(di ) , where the sum is extended to the four h - bonds . the globally funneled energy landscape theorized for protein folding allows for a significant reduction in the complexity of the description provided that simple reaction coordinates are found to capture the folding process.(1 ) in this study , the simulation trajectories were projected along the observables hb and rmsd and binned into n two - dimensional microstates . thus , each bin / microstate is defined by both a range of hb values and a range of rmsd values . in this way , a discretized trajectory s(t ) n]t  [ 0,tf ] was obtained , where tf is the total simulation time . uniform binning ( i.e. , all bins having the same width along each of the two observables ) resulted in poor statistics for some of the transitions between microstates , which compromised the accuracy of the rate matrix diagonalization ( data not shown ) . thus , a variable binning strategy was adopted , which requires a few iterations of the following algorithm . starting from a single two - dimensional bin containing all the data points , bins containing more than ntb = 12  000 points were evenly split into four smaller identical bins with halved width along each dimension , and the data points were redistributed in the four smaller bins . no splitting was performed for bins whose width was equal or smaller than ( hb = 0.25  rmsd = 0.625 ) , both to limit the total number of bins and to prevent artificial splitting of very similar conformations . this procedure favors a more uniform distribution of the data points into the bins ( i.e. , the microstates ) , and eventually leads to more accurate estimates for transition probabilities between microstates . the binning procedure was applied separately for the simulations at 300 , 330 , and 360 k , resulting in different numbers of bins / microstates at each temperature . the observed probability pi of each microstate is given bywhere t is the time and  is the kronecker delta function . the number of transitions observed from microstate i to microstate jnji after a lag time t is thus expressed as below , the msm is derived following the formalism proposed by hummer and co - workers.(18 ) assuming that , after the initial lag time t , the dynamics observed along the simulation trajectories is markovian , the probability pi satisfies a master equation : where the vector p(t ) represents the probability for each microstate at time t and r is a transition rate matrix with constant elements rij  0 for i  j , rii  0 , and irij = 0 . equation 3 has the solution if the system converges to an equilibrium distribution , the equilibrium probability p is a right - eigenvector of r with zero eigenvalue , and the system must satisfy detailed balance so thatthe matrix exponential of eq 4 can be calculated through the diagonalization of the symmetric matrix r:where p is the diagonal matrix with elements ( pieq ) . for the diagonalization , the orthogonal matrix u of eigenvectors of r must satisfy the relation ru = u , where  is the diagonal matrix with eigenvalues of i . since uu = uu , where u is the transpose of u , eqs 4 and 6 can be combined : the left and right eigenvectors of r can be obtained from the eigenvectors of r ( note that they have the same eigenvalues ) : in addition to the stationary solution p with zero eigenvalue , the other right - eigenvectors of r have negative eigenvalues representing decay rates , and the vectors elements sum up to zero . so a system in state i at time t = t0 has a probability p(j , t0 + t|i , t0 ) of being in state j at time t0 + t given by the rate matrix r was estimated using bayesian inference by constructing a posterior distribution of the model parameters from the simulation data . if a uniform prior distribution of the model parameters is assumed , the posterior distribution is proportional to a likelihood function l , which gives the probability of observing the same motions on the hb - rmsd reaction coordinates in the markov model as those seen in the simulations : this may be rewritten in the form a monte carlo sampling of the values of r and p compatible with the simulation data was obtained using ln(l ) as the energy function and by constructing the initial rate matrix r from the equilibrium probabilities pi ( eq 1 ) and the observed transitions nij ( eq 2 ) after symmetrization and normalization . a simulated annealing strategy was used where the monte carlo temperature , initialized at 10 , was gradually lowered to 1 following an exponential decay as a function of monte carlo steps . only the data obtained with monte carlo temperatures below 1.005 were retained for the sampling . the monte carlo step size was dynamically adapted in order to obtain an acceptance ratio between 10 and 30% . a total of 10 monte carlo steps were performed , and samples of r and p values were saved every 1000 steps . each sampled r and p pair was individually used to compute the eigenvalues and eigenvectors of the associated markov model . eigenvalues and eigenvectors of r were averaged over the sampled data . on the basis of perron clustering theory,(31 ) the microstates were then partitioned into macrostates using information from the slowest non - null left eigenvectors of r. because of the gap between the third and fourth non - null eigenvalues , we considered the three slowest eigenvectors for partitioning . each macrostate is defined as the set of microstates having the same signs of the components of the three slowest left eigenvectors used for partitioning . the leader algorithm(33 ) was used to group structures whose distance root - mean - square deviation between c and c atoms was less than 0.8  . this value for the clustering cutoff was chosen as a compromise between smaller values , yielding many small clusters and poor intercluster statistics , and larger values , resulting in heterogeneous structures being grouped in the same cluster . the clustering and other structural analyses were performed using the program wordom.(34 ) the commitment probability pcommit, for a protein conformation is the probability that it will fold into a given state  before reaching state . there are several methods to measure the pcommit of a conformation.(35 ) one method is to start a large number of short md simulations from the given conformation , and record the number of trajectories that commit to a given state of the molecule.(36 ) the number of trial simulations required for a reliable evaluation of pcommit for the many conformations saved along a trajectory makes this method unfeasible . an alternative method estimates commitment probabilities for all structures in an equilibrium folding / unfolding trajectory without the need for any additional simulations . the probability pcommit, of a single structure i reaching either of the states  or  is given bywhere n is the number of md simulations starting from structure i which reach state  and ntot is the total number of committed md simulations starting from i. since structurally similar conformations have similar kinetic behavior , they also have similar pcommit,.(11 ) structures along a trajectory were grouped in structurally similar clusters , as described above , and the segment of md trajectory following each structure was checked for reaching either the states  or . the pcommit, of a cluster c is defined as the ratio between structures folding into state  and the total number of committed structures in the cluster . each structure in the cluster although the msm identified several slow relaxation processes in tz1 , we are mostly interested in the one between the folded and the denatured region of the peptide , i.e. , the main folding transition , which , in the present case , is represented by the slowest non - null eigenmode of the rate matrix . thus , we considered only the sign structure of the slowest non - null eigenvector and we selected the three most populated microstates with a positive eigenvector component ( i.e. , belonging to the folded region ) and the most populated microstate with a negative eigenvector component ( i.e. , belonging to the denatured region , the denatured region has a much larger population than the native region so a single microstate is sufficient ) and used them as targets  and  for measuring pcommit . the  and subsets are meant to represent the bottom of the free energy well of the native and denatured region , respectively . the ability of the approximate commitment probability as a reaction coordinate to locate the dominant transition between the folded and denatured states can be evaluated by examining the conformations on a transition pathway . following hummer,(37 ) a transition pathway ( tp ) is defined as any section of a trajectory leaving a state a and proceeding directly to state b ( i.e. , without crossing back into a ) , and vice versa . in our calculations , a and b are the conformations with pcommit  0.2 and pcommit  0.8 , respectively . a minimum path length of 0.08 ns was imposed to rule out rapid fluctuations at the boundaries of the tps . the conditional probability that a conformation with a given value of the pcommit coordinate will be on a tp , p(tpab|pcommit ) , is then measured directly from the trajectory , and its distribution is compared with the limiting case of the pure unidimensional diffusion to ascertain the quality of our approximate pcommit as a reaction variable for the folding transition . the conformation space of the peptide , as sampled by molecular dynamics , can be mapped to a network describing the significant free energy minima and their dynamic connectivity , without requiring arbitrarily chosen reaction coordinates . some studies(6 ) have found poorly chosen reaction coordinates can mask important details , such as heterogeneity of denatured or transition states , especially in the case of potentials based on physiochemical principles , like charmm . in this case , the  nodes  of the network are the clusters defined previously by the cluster analysis , with direct transitions between clusters along the trajectory representing the  edges  of the network.(9 ) two clusters are connected by an undirected link if either they include a pair of structures that are visited within 20 ps or they include structures that are separated by one or more conformations belonging to clusters with less than 20 structures each . the nodes of the network were ordered on the graph so that close nodes correspond to clusters with a similar link pattern , and colored according to the markov macrostates defined above . the stable states of the peptide were identified by optimizing a msm using bayesian analysis , as described in initial calculations carried out using only the backbone h - bonds for the definition of the microstates did not show markovianity , as reported by the eigenvalue spectrum which did not converge at any lag time ( data not shown ) . simulations were then projected along both hb and rmsd , and the data were binned into a discrete number of two - dimensional microstates defined by both a range of hb and rmsd values ( figure 2a ) . thus , 55 , 53 , and 41 microstates were obtained at 300 , 330 , and 360 k , respectively . the variable width of the microstates / bins , obtained using the adaptive binning procedure described in the  methods  section , reduces the presence of scarcely populated bins which would introduce large uncertainties on the elements of the transition matrix . the ensemble of rate matrices compatible with the simulation data was obtained using the bayesian analysis and the likelihood function l ( eq 11 ) . the rate matrices were then diagonalized . at all three temperatures ( 300 , 330 , 360 k ) , the slowest non - null eigenvalue converged for lag times not less than 2.56 ns , and a well - defined gap was observed between the slowest and second - slowest non - null eigenvalues ( the eigenvector corresponding to the null eigenvalue is the equilibrium state ) and between the third and fourth non - null eigenvalues ( figure 2b ) . the relaxation times for the three slowest processes at the three temperatures are reported in table 1 . although the slowest relaxation process at 360 k is on average slightly faster than the lag time , the subsequent analysis of the corresponding macrostates provided results in line with those obtained at the lower temperatures , so we decided to present the 360 k data as well . ( a ) equilibrium distribution of the microstates of the markov state model fitting simulation data at 300 k ( left ) , 330 k ( middle ) , and 360 k ( right ) . the variable width of the microstates in the hb / rmsd space is obtained via an adaptive binning algorithm ( see  methods  ) to limit the presence of poorly populated microstates . ( b ) relaxation times ( inverse eigenvalues ) of the markov state model as a function of lag time for the simulations at 330 k ( error bars at 99% confidence interval ) . first ( left ) , second ( middle ) , and third ( right ) average non - null left eigenvectors of the rate matrices at ( c ) 300 k , ( d ) 330 k , and ( e ) 360 k. the color scale ranges from black for the smallest ( most negative ) components to white for the largest ( most positive ) components , with the midpoint of the scale fixed at 50% gray for the zero components . the sign of the components of an eigenvector of the rate matrix ( note that components of left and right eigenvectors of r have the same sign ) tells us which microstates exchange population during the corresponding relaxation process ( e.g. , the microstates with positive components of the eigenvector lose population in favor of those with negative components , or vice versa , depending on the initial conditions ) . on the basis of this idea , which has been rigorously formalized in perron clustering theory,(31 ) the signs of the components of the first three non - null left eigenvectors of the rate matrices , which give information on the slowest , second - slowest , and third - slowest relaxation processes , respectively , were used to identify four significantly populated stable macrostates , as detailed in inspection of the elements of the first eigenvector reveals that it describes the relaxation process between the folded and denatured regions of the peptide s conformational space . at 300 k ( figure 2c , left ) , the folded region is roughly defined by a maximum rmsd of 6  and between 0.5 and 4 formed hb contacts , and the denatured region is defined by a maximum of 0.25 formed hb contacts for 2.25  < rmsd < 6.25  and a maximum of 1 formed hb contact for rmsd 's above 6.25  . the second eigenvector yields more detailed information about the denatured region , which at 300 k ( figure 2c , center ) is comprised of two distinct macrostates occupying separate regions on the hb - rmsd projection : the first macrostate is roughly defined by hb < 0.25 and 2.25  < rmsd < 3.75  and will be referred to as state d1 . the second macrostate is roughly defined by hb < 1 and rmsd > 3.75  and will be referred to as state d2 ( at 360 k , in figure 2e , middle , there is some scatter with a microstate from d1(dark ) in the d2 region ; however , this microstate has an extremely low population as reported by figure 2a , right , and its presence does not affect the validity of our statement ) . the most native - like subset will be referred to as the native state n and will comprise roughly conformations with hb > 0.5 and rmsd < 2.25  . the intermediate region ( at 300 k with roughly 0.25 < hb < 2 and 2.25  < rmsd < 6.25  , figure 2c , right ) will be referred to as state i. the eigenvectors at 330 and 360 k show a very similar picture , as illustrated in figure 2d and e. in figure 3 , information from all three eigenvectors is merged to yield a picture comprising four stable macrostates at all three temperatures : n ( magenta ) , i ( yellow ) , d1 ( blue ) , and d2 ( green ) . markov macrostates at ( a ) 300 k , ( b ) 330 k , and ( c ) 360 k obtained by grouping together microstates with the same signs of the components of the three slowest eigenvectors of the msm . the native state n is shown in magenta , the intermediate state i in yellow , the denatured state d1 in dark blue , and the denatured state d2 in green . conformational clusters correspond to the nodes of the network , and the observed transitions correspond to the links between nodes . this mapping allows for the definition of a vector for each cluster whose components are the network distances ( i.e. , the minimal number of separating links ) between the cluster and a set of reference clusters ( e.g. , the most populated ones ) . the network can be projected onto a two - dimensional graph , using the first two principal components of the covariance matrix of the cluster vectors.(38 ) this representation has the advantage of preserving the kinetic similarities between clusters . thus , the network could be represented in a plot where the clusters were colored according to their most populated state , identified by the msm analysis , with the same color scheme as in figure 3 . the network analysis performed on the simulations of tz1 at each temperature ( see figure 4 ) show the n and i states occupying a small , localized region of the network while the denatured states d1 and d2 occupy the majority of the network . at all three temperatures , the mapping of the macrostates from the msm analysis onto the kinetic network shows that the markov states are kinetically partitioned , i.e. , different states occupy different regions of the network , with limited overlap . kinetic networks describing the conformation space of the tz1 peptide at ( a ) 300 k , ( b ) 330 k , and ( c ) 360 k. the nodes of the network are colored according to their most populated macrostate , with the n state being magenta , the i state being yellow , the d2 state being green , and the d1 state being blue . the network pictures are obtained using the program visone.(38 ) the localization of the n state is particularly prominent at 300 k , where the relative depth of the free energy minima is larger with respect to temperature , as reported by the low link density region separating the native clusters from the rest of the network in figure 4a . this confirms the presence of free energy barriers between the markov macrostates , and that the observed denatured states are in a different free energy basin than the folded states . in agreement with the kinetic measurements , the presence of the gap is less well - defined at 330 and 360 k , corresponding to a smaller free energy barrier relative to kbt , separating n from i. we note that the increase in link density between separate free energy basins as the temperature increases is equivalent to an increase in the available pathways used by the system to connect the basins , in analogy to the findings for the trp - cage mini protein.(39 ) at all three temperatures , there is also evidence of a free energy barrier between the d1 state and the most populated region of the d2 state ( green ) , although regions where the two states are mixed are also present on the network . each relaxation process can be associated with the crossing of a free energy barrier defined on the free energy landscape of the peptide . the highest barrier , corresponding to the slowest relaxation process , separates variably folded microstates from denatured microstates . the folded microstates , which share the same sign of the first eigenvector components ( light - colored microstates in figure 2c  e , left ) , comprise both microstates in n and i. the denatured microstates ( sharing the opposite sign to the folded microstates for the elements of the first eigenvector , dark microstates in figure 2c  e , left ) comprise both d1 and d2 . smaller free energy barriers split the folded microstates in i and n and the denatured microstates in d1 and d2 . disregarding the smaller free energy barriers , we specifically looked for the main folding transition state associated with the highest free energy barrier . approximate commitment probabilities for the main folding transition were measured for all the structures collected along the trajectories at each temperature , using the most populated microstates pi in the folded and denatured regions ( as determined by the sign of the first non - null eigenvector ) as target states ( see  methods  ) . the distribution of commitment probabilities ( figure 5a ) shows two distinct peaks at pcommit  0.0 and pcommit  1.0 , which correspond to the denatured and folded states , respectively ( the distribution of pcommit for intermediate values is noisy due to the relatively low number of conformations populating this region ) . the approximate procedure used to measure the commitment probability , which depends on the clustering algorithm and sampling , may in principle compromise its reliability as a reaction coordinate . hummer and collaborators proposed a method to evaluate the quality of a reaction variable r by measuring the conditional probability , p(tp|r ) , that a conformation will be on a transition pathway ( tp ) for a given value of r. he showed(37 ) that , in the limit of diffusive dynamics along r , the transition path probability p(tp|r ) is related to the exact commitment probability (r ) by p(tp|r ) = 2(r)(1  (r ) ) . thus , p(tp|r ) in the limiting case of a pure diffusive system has a maximum of 0.5 at the transition state , i.e. , for (r ) = 0.5 . this finding has been used to assess the quality of reaction variables for describing the dynamics of a generic system . in those studies , the authors compared the p(tp|r ) obtained using several reaction variables . those variables which provided p(tp|r ) distributions close to the diffusive limit ( i.e. , peak close to 0.5 , bell shaped distribution around the peak ) are considered good , as they can separate stable and reactive species . vice versa , those variables giving rise to broad p(tp|r ) distributions with a small peak ( < 0.3 ) can not be used to tell unequivocally if a structure is reactive ( i.e. , if it belongs to the transition state ) ; thus , they do not act as good reaction coordinates . in the present case , the distribution of p(tp|pcommit ) , shown in figure 5b , is bell - shaped with a maximum of 0.48 at 300 k , 0.49 at 330 k , and 0.48 at 360 k centered on the region 0.4  pcommit  0.55 , almost reaching the ideal limit of 0.5 observed in the case of pure unidimensional diffusion . this demonstrates that the approximate pcommit reported in the present work represents a good reaction coordinate . ( a ) distribution of approximate commitment probabilities shown at 300 k ( solid line ) , 330 k ( dashed line ) , and 360 k ( dash - dotted line ) . ( b ) the conditional probability of being on a transition pathway for a given value of the approximate commitment probability is shown at 300 k ( solid line ) , 330 k ( dashed line ) , and 360 k ( dash - dotted line ) . once the metastable macrostates and transition state were identified by the bayesian and commitment probability analyses of the markov model , respectively , their structural and thermodynamic properties could be examined . the populations of states d1 , d2 , i , n , and t were used to calculate their relative stabilities ( table 2):where p(a ) and p(b ) are the populations of the two states , na is the avogadro constant , and kb is the boltzmann constant . the average internal potential energy of each macrostate epot was calculated from the internal potential energies recorded for each structure along a trajectory ( table 3 ) . the relative enthalpies of the denatured and native state were also compared : where h(t)u n is the enthalpy change at temperature t and epot(u) and epot(n) are the average internal potential energies of the denatured ( d1 + d2 ) and native ( n ) state , respectively . the average number of backbone h - bonds and average rmsd 's for each state are shown in tables 4 and 5 . as shown in table 3 , the n state has the lowest average potential energy and the denatured states d1 and d2 have the highest average potential energies at 300 , 330 , and 360 k , as expected . at all three temperatures , the enthalpy increases in the denatured state relative to the native state ( h(t)u  n = 4.1  4.9 , 4.7  5.5 , and 4.2  5.9 kcal / mol at 300 , 330 , and 360 k , respectively ) , although the error on the difference is quite large . a structural comparison of the various states of the peptide the secondary structure of the peptide according to stride(42 ) was computed along the trajectories using vmd(43 ) for the macrostates d1 , d2 , i , n , and t. for the compact denatured state d1 ( figure 6 , column 1 ) , the hairpin is quite  open  , and with the h - bonds between residues e5 and k8 only loosely formed ( hb < 0.01 ) , the shape of the backbone fluctuates significantly within this state . the turn sequence extends from residues w4 to n7 at the three simulation temperatures , while the rest of the residues are in the coil region . the four tryptophan side - chains appear to stabilize this state . in this case , the tryptophan residues appear to minimize their exposed surface area by packing between the backbones of the two strands . the second denatured state d2 ( figure 6 , column 2 ) has , on average , a greatly extended turn sequence from residues w2 to t10 at 300 , 330 , and 360 k , though the secondary structure is less homogeneous than for the other macrostates . the hairpin is more open with much higher fluctuation in both backbone shape and the positions of the tryptophan side - chains . the i state ( figure 6 , column 3 ) , in which the native g6n7 turn sequence is preserved , features -bridge structure at residues e5 and k8 and additional turn structure for residue w4 at 300 , 330 , and 360 k. the rest of the peptide is random coil . the hairpin remains relatively  closed  , and the tryptophan residues form a disordered hydrophobic core on one side of the hairpin . the shape of the backbone shows slightly less variation within this macrostate , primarily toward the terminal residues , since the turn sequence is stabilized by a strong h - bond between residues e5 and k8 . the n state conformations ( figure 6 , column 4 ) at all three temperatures are characterized by the preservation of the g6n7 turn sequence and the presence of extended -sheet structure in the four residues on either side of the turn , residues w2e5 and k8w11 . stabilized by the presence of five strongly formed backbone h - bonds , the ensemble of conformations shows very little variation in the shape of the backbone within the n state , and the tryptophan residues are regularly packed in the center of the hairpin . conformations of the four stable macrostates identified by the markov state model analysis ( d1 , d2 , i , and n ) and the transition state identified by the commitment probability ( t ) at ( a ) 300 k , ( b ) 330 k , and ( c ) 360 k. the unstructured random coil is shown in gray , turn sequence in tan , extended -sheet structure in red , and -bridge structure in blue . several backbone conformations for each state are superimposed to illustrate the conformational variability within the state , and representative positions of the four tryptophan side chains are shown . the c atoms of the conformations were aligned using the first nmr conformation as a reference . at 300 and 330 k , the transition state t ( figure 6 , column 5 ) identified from the commitment probability in figure 5 is structurally very similar to the intermediate state , though the variation in structures is greater for this less energetically stable state . at 360 k , while the g6n7 turn sequence is unchanged , the -bridge structure of residues e5 and k8 dominant at the two lower temperatures is replaced by two residues of the extended -sheet ( w4e5 and k8w9 ) . the fluctuation in the backbone shape is also similar to that of the intermediate state . at 300 k , there are two formed h - bonds between residues e5 and k8 , while at the two higher temperatures only the h - bond closest to the turn is formed . the changes in the secondary structure content between the states along the peptide sequence were also examined . within each of the d1 , d2 , i , n , and t states , the numbers of conformations with random coil , turn , extended -sheet , or -bridge structures were calculated as a fraction of the total number of conformations for each of the 12 residues at the three simulation temperatures . since the changes in secondary structure between the states at each temperature are nearly identical ( see figure 6 ) , only the 330 k data are shown in figure 7 . residue - specific secondary structure distribution for the states d1 ( blue ) , d2 ( green ) , i ( yellow ) , n ( red ) , and t ( black ) at 330 k for ( a ) the turn sequence , ( b ) extended -sheet , ( c ) -bridge structures , and ( d ) random coil . these data highlight the homogeneity of the n state , with all the conformations having a similar secondary structure , as expected . while the d1 state is also structurally homogeneous ( with the exception of residues k8w9 ) , the outer tryptophan pair ( w2 and w11 ) and the two threonine residues ( t3 and t10 ) in the d2 state can variably be in either a turn conformation or a random coil . similarly , the intermediate and transition states are less structurally homogeneous , particularly for the residues on either side of the g6n7 turn sequence ( w4e5 and k8w9 ) , for which no secondary structural feature is present in more than half of the conformations within the state . it is also noteworthy that in d1 not only g6 and n7 show recurrent turn secondary structure , but also w4 and e5 ( figure 7a ) . in addition , -sheet formation in d1 , although relatively rare , involves residues w2w4 and residues e7w9 ( figure 7b ) . these two pieces of data suggest the presence of out - of - register hairpins in d1 with the turn shifted one residue to the n terminus . the relevance of all formed contacts was determined by calculating structural -values for each residue in the intermediate and transition states at the three simulation temperatures . in each individual structure , a van der waals contact was defined when the distance between two atoms was less than 6.0  , and was defined as  formed  if it occurred for at least 2/3 of the structures within a state . nearest - neighbor contacts ( both interatomic and between residues ) and backbone atoms were excluded from the calculation . if the fraction of n - state structures in which the contact i is formed is n(i ) and the fraction of structures in a state x in which the contact i is formed is x(i ) , the structural -value for state x iswhere nc(r ) is the group of contacts occurring in conformations of the native state for each residue r. snat(r ) usually takes values ranging from 0 to 1 ( although values out of this range are also possible ) , 0 indicating that the structure around residue r is conformationally far from the n state and 1 indicating that the structure is conformationally similar to the n state . as can be seen in figure 8 , the central residues of the peptide are most native - like , with a general decrease in nativeness toward the terminal residues . the number of contacts formed in the intermediate state is largely independent of temperature , while the transition state becomes less native - like as the temperature is increased . at the lower temperatures , the t state snat values are larger than those in the i state , indicating a more native - like packing of the side - chains . figure 8 confirms that the g6n7 turn sequence is strongly native - like for both the intermediate and transition states at all three temperatures . structural -values per residue in the i and t states at 300 , 330 , and 360 k. the optimized markov state approach employed here is based on several assumptions whose validity must be considered . the first of these is that the dynamics of the system is markovian , or memory - less . this implies that the future dynamical evolution of the conformation of a molecule depends only on its current state and not on the previous history of the molecule . in practice , the dynamics of a system is unlikely to be completely markovian at short time scales,(19 ) so the mean lifetime of the stable states of a system must be long compared to the length of the memory , allowing fast initial non - markovian dynamics to be neglected . a preliminary test of the markovian assumption is based on the convergence of the spectrum of the eigenvalues of the rate matrix , as the lag time is varied . this test showed us that a representation of the system obtained by binning the trajectory along a single observable ( the number of native backbone hbs ) did not provide markovian dynamics ( as it did in other instances ) . this means that the information conveyed by the hb is not sufficient to describe the kinetics of folding of the peptide . the method was therefore extended to use both hb and rmsd ( see methods ) , and in this case , the slowest eigenvalues at the three simulation temperatures studied were found to converge for lag times larger than 2.5 ns ( figure 2b ) . as a further and more stringent test of the markovian assumption , the decay times of the autocorrelation functions of the first three eigenvectors of r projected along the simulations were examined . figure 9 shows that these decay times are in good agreement with the decay times expected from the corresponding eigenvalues ( table 1 ) , confirming that this representation is markovian . the markovianity of the representation provides a validation , a posteriori , for both the choice of observables ( hb and rmsd ) and the fine graining of the binning procedure . indeed , if either of the two was not good enough , then the representation would not be markovian . we note also that other choices of observables should provide similar slow eigenmodes , as long as the observables report the same amount of information about the slow dynamics of the system.(17 ) autocorrelation functions for the first three eigenvectors ( continuous lines ) of the rate matrix r at 330 k and the corresponding exponential decay associated with the eigenvector ( dashed lines ) . the eigenvector analysis identified three significant relaxation processes , the slowest being the one between the folded and denatured states of the peptide , with relaxation times ranging from 2 to 20 ns , according to the temperature ( table 1 ) . the relaxation times measured in this study can not be directly compared with experimental data , as the simulations were deliberately run in the absence of solvent viscosity to maximize the number of sampled transition events . however , the 2 orders of magnitude reduction in the transition times measured at 300 k ( 20.0  0.6 ns , table 1 ) with respect to the experimental data ( kobs1 = [ kf + ku ] = 4.7  0.3 s at 296 k(36 ) ) is in line with the reduction observed in a variety of other systems ranging from helices to three - stranded -hairpin peptides simulated in the same way , indicating that the force field yields a separation of time scales consistent with experimental data . interestingly , the second slowest relaxation , with relaxation times ranging from 1.4 to 4.6 ns ( table 1 ) , occurs between a highly populated , compact state ( d1 ) and a less densely populated state encompassing a broader range of hb and rmsd values ( d2 ) . evidence of this separation of the d1 and d2 states was also indicated in the kinetic network analyses by regions of low link density separating the two ( figure 4 ) . the distance from a pure exponential relaxation for the autocorrelation function of the associated eigenvector ( figure 9 ) and the presence of a limited overlap ( figure 4 ) in the kinetic network projection between d1 and d2 states , as reported by the presence of network regions containing clusters from both states , indicates , however , that even the rmsd and hb may be flanked by some other observable to improve the separation between the d1 and d2 states . we also note that , because d1 and d2 may not have been perfectly separated by the hb - rmsd - based msm , the relaxation times between the two states may be underestimated . these relaxation times are already close to those measured for the main folding transition ; thus , we can not completely exclude the possibility that the d1/d2 exchange occurs via the native state , or , in other words , that the native state represents a kinetic hub connecting the two heterogeneous denatured states as recently proposed for 3s(12 ) and for villin and ntl9.(49 ) on the other hand , here the native state is not very stable at any of the simulated temperatures ; thus , our system is in a regime where the kinetic hub picture may not necessarily hold . it is plausible that the lack of markovianity of the representation based only on hb is due to the complexity of the denatured state and the presence of two distinct free energy basins d1 and d2 , which can not be separated by hb alone . the third significant free energy barrier at 330 and 360 k was found between the n and i states ( figure 2d , e , right ) , with relaxation times ranging from 0.8 to 1.8 ns ( table 1 ) . this barrier , which is visually more evident from the hb - rmsd projections than the one between d1 and d2 , is actually smaller than the latter , as it allows for a faster relaxation ( i.e. , larger eigenvalue ) . this is a clear demonstration of the need for a careful kinetic analysis to correctly discriminate the populated macrostates of a biomolecular system , especially when the nature of the stable states is unknown a priori , as in the present case , or in any attempt to characterize the structures of the denatured state . at 300 k , the third relaxation process primarily involves a population exchange between i and d2 , with only a weak participation of n ( figure 2c , right ) , in contrast to what happens at higher temperatures . we note , however , that this result is still sufficient to separate n from i. the disadvantage of using a msm based on projections of the trajectory upon simple observables ( i.e. , hb and rmsd ) is that the kinetic information contained in the simulation data may be partly lost , an issue that becomes particularly serious when kinetically unstable transition states are considered.(6 ) the decomposition in macrostates provided by the msm does not allow for the direct identification of reactive species . indeed , the distribution of the probability of folding for conformations belonging to the i state shows that i is enriched in reactive conformations with respect to the overall trajectory . however , it also shows that the largest majority of the conformations in i are not reactive ( figure 10 ) . then , the use of the kinetic network analysis , which allows for the definition of a reliable reaction coordinate , pcommit , becomes necessary for the correct identification of the reactive species . probability density of the commitment probability distribution for state i ( continuous line ) and for the overall trajectory ( dotted line ) . the denatured state of tz1 in the present simulations appears heterogeneous , because the two substates d1 and d2 have a relatively slow interconversion rate . the state d1 is a compact state where the turn at residues g6 and n7 is formed , while the two strands of the hairpin are not paired by backbone hb and interlocking trp side chains like the native state but by a disordered arrangement of the trp side chains in a conformation that possibly minimizes the nonpolar exposure to the solvent . the state d2 , on the other hand , is more diffuse and heterogeneous with a lower amount of residual structure , although it retains a certain tendency to form a turn close to residues g6 and n7 ( figure 7 ) . it is worth pointing out that the compact structures with buried tryptophan side chains observed here for the states d1 , and to a lesser extent d2 , may be the result of an overstabilization of hydrophobic interactions due to the approximate model for the solvent . this may also explain the low population of the native state observed here even at 300 k , in contrast with the experimental midpoint transition temperature of 323 k.(25 ) although the comparison with explicit water simulations of the same peptide(22 ) is made difficult by the different temperatures used in that study , no strong propensity to sample compact conformations in the denatured state was observed in those simulations . nor were they able to distinguish the two components d1 and d2 of the denatured state , although a further investigation ( g.s . unpublished data ) revealed that the tendency to form a turn in position 6 and 7 in the denatured state was also present in the explicit water simulations . high temperature explicit water simulations of tz2,(50 ) a peptide similar to tz1 but with a different turn sequence , and an analysis using an automatic markov state decomposition(13 ) revealed considerable structure in the denatured state in analogy with what has been found here , although in those studies structure in the denatured state appeared to be stabilized by non - native backbone h - bonds while here it is predominantly stabilized by hydrophobic interaction between tryptophan side chains . those studies also pointed out the presence of out - of - register hairpin states , some of which are compatible with the n - terminal shift of the turn observed here in some of the conformations of d1 . the presence of the significantly populated intermediate state ( i ) in the current simulations was also observed in explicit - solvent high temperature molecular dynamics simulations of the same peptide,(22 ) although the detailed structures differ . indeed , in the present simulations , state i is stabilized by a hydrogen bond between k8 and e5 , while , in the explicit water simulations , a hydrogen bond between t3 and t10 was also formed . previous simulations performed on the related tz2 peptide led to contrasting results . while large numbers of short implicit - solvent md simulations suggested a two - state system,(36 ) other simulations of tz2 using implicit solvent models(51 ) favored a noncooperative folding hypothesis , even if the projection of the free energy landscape showed a minimum with characteristics similar to the i state identified in this study . explicit solvent simulations of tz2 also revealed the presence of an intermediate state similar to the i state . experimental characterization of tz2 by isotope - edited 2d ir spectroscopy(55 ) revealed the presence of partially folded states with only the midrange native h - bonds formed along the hairpin , a bulged turn , and variably frayed termini . these conformations resemble the i state identified for tz1 here and , to a larger extent , that observed in the explicit water simulations.(22 ) explicit water simulations of the hairpin b1 from protein g,(56 ) which has a similar sequence to the trpzip peptides , reported the presence of at least one kinetic intermediate with characteristics similar to the one observed here . the presence of an intermediate state fits the available experimental data on the folding of tz2 and related peptides , where the unfolding rates are mainly determined by the robustness of the hydrophobic core packing and the folding rates are determined by the turn propensity of the sequence.(57 ) the transition state t of tz1 is characterized by a well formed turn region and the proximal part of the strands in extended -sheet conformation ( figures 6 and 7 ) . with respect to secondary structure content , t is similar to the i state although slightly less native - like . in terms of side chain packing , the t state presents a compact and disordered arrangement of the tryptophan side chains similar to i , although the number of native side - chain contacts in t is larger than in i ( figure 8) , indicating the importance of correct side - chain packing for the folding transition to complete . this picture of the transition state is similar to the one obtained from explicit solvent simulations(22 ) where the rate limiting step for folding was identified in the correct packing of the tryptophan side - chains , while proximal h - bonds and the turn region were already formed at transition . the structure of the folding transition state was found to be relatively insensitive to changes in temperature in proteins like ci2 and protein a.(60 ) however , smaller structural units such as those represented by small peptides may be more sensitive to environmental changes . while explicit water simulations of the unfolding of the tz1 peptide(22 ) showed subtle structural changes in the transition state structures as a function of temperature , the high temperature and lack of proper thermodynamic equilibrium in the simulations prevented a deeper exploration of these phenomena . here , we observe that the similarity of the t state to the native state n , as reported by the structural  values and by the superposition of conformations of t in figure 6 , decreases with temperature . this , however , should not be interpreted as anti - hammond behavior because this comparison does not take into account the corresponding changes in the denatured state . indeed , the conformational distribution of n and d1/d2 is also temperature dependent . experimentally , movements of transition states relative to the ground states are measured in terms of the tanford beta value,(61 ) which is thought to approximate the ratio between the solvent accessible surface area buried in t and in n , upon folding from d. accessible surface area calculations on the conformations identified here as belonging to t , n , and d2 ( table 6 ) show small differences between the states , relative to the spread of the distributions . this leads to large errors on the determination of the tanford beta value , which rules out drawing a definitive conclusion about temperature dependent transition state movements . a possible reason for the large error lies in the fact that d2 in our simulations , although less compact than d1 , is still a relatively compact state ; thus , the change in solvent accessible surface area upon folding is limited . here , we showed how the combined use of msm analysis and kinetic network analysis allows for a detailed , concise , and consistent characterization of both the metastable and the transition states for folding of peptides . besides states identified in earlier simulations , the markov state analysis of trajectories projected on standard observables allowed for the identification of different subspecies in the denatured state of the tz1 peptide , which could otherwise be easily missed by standard inspection of the projected trajectory . the metastable states as identified by the markov state analysis provided a solid base to perform approximate commitment probability calculations using the kinetic network representation of the peptide dynamics and correctly identify the reactive species . the approximate commitment probabilities were proven to represent a good reaction variable . the markov state analysis and the kinetic network analysis also provided a consistent reciprocal validation of the results , showing that the markov macrostates are kinetically localized and well partitioned along the network . the picture of the folding process of tz1 emerging from the analysis is the following . the slowest conformational change observed in tz1 is the conversion between folded and denatured species . the other slow processes involve exchange between a compact and a less compact subspecies of the denatured state and exchange between the fully native and partially folded intermediate state with an h - bond formed close to the turn . the transition state for the slowest relaxation process ( i.e. , folding ) is characterized by the formation of h - bonds proximal to the turn and a well structured turn region , while packing of the tryptophan side - chains is not complete , although it encompasses more native contacts than in the intermediate state .