entropy is a key property to understand a wide variety of physical , chemical , and biochemical phenomena .
however , the estimation of absolute entropies and entropy differences from computer simulations is a long - standing problem and one of the current challenges in computational chemistry .
the calculation of reliable absolute entropies from molecular dynamics ( md ) simulations is intrinsically difficult because the absolute entropy is a measure of the overall extent of phase space ( ps ) accessible to a molecular system .
however , absolute single - molecule entropies can be estimated based on an analytical approximation to the configurational probability distribution corresponding to the ps accessed by a simulated system.(2 ) the underlying theory , assumptions , approximations , and alternative practical implementations have been recently reviewed .
the relationship among quasi - harmonic ( qh ) , essential - mode , and normal - mode analyses has also been investigated.(11 ) for an extensive review of the subject , not limited to the qh approach , see also refs ( 11 ) and ( 1619 ) and references therein . the difference between the true entropy of a simulated system and
its qh estimate arises from ( i ) anharmonicities ( i.e. , non - gaussian behavior ) in the probability distributions along individual qh modes and ( ii ) correlations among the probability distributions associated with different qh modes ( beyond the pairwise linear correlations accounted for ) .
these effects are neglected in standard qh analysis and ( nearly ) always lead to a negative entropy contribution.(11 ) a method to correct for both artifacts was recently described.(11 ) point ii is of particular relevance when trying to estimate entropy differences between two conformational states of a molecular system because error cancellation can not be guaranteed a priori . by taking into account correlation effects of increasing order , entropy estimates based on corrected qh analysis aim at capturing the entropy corresponding to the entire ps sampled ( see figure 1 in ref ( 11 ) ) .
thus , this approach is not limited to systems sampling single qh basins and allows capturing conformational transitions .
the qh entropy upper - bound sqm , o ( eq 10 with eq 12 ; empty bars ) and the improved entropy estimate s ( eq 20 ; hatched areas ) are displayed scaled by the number n of system particles .
the relative ( % ) values of the cumulative correction term fclcrc ( eq 21 ; bar labels ) are reported as a measure of the importance of the deviation from the qh approximation . from left to right : two -peptides in methanol at high temperature ( f , folded ; u , unfolded ; a , all ; t , 298 instead of 340 k),(11 ) the 11 glucose - based disaccharides ( w , free in water),(28 ) dipalmitoylphosphatidylcholine , dppc ( i , inserted in a hydrated bilayer),(24 ) the cavity and its gating loop of the w191 g mutant of cytochrome c peroxidase ( k , bound to a k ion with closed gating loop ; b , bound to 2-amino-5-methylthiazole with closed gating loop ; o , bound to a k ion with open gating loop),(20 ) the h - ras lipopetide anchor,(27 ) and the cc-peptide ( this study ) . for the disaccharides,(28 ) corresponding mean entropy values are displayed ( a vertical bar represents the range of values ) .
see supporting information for details . in the present article , we expand the previous study in ref ( 11 ) . a general formulation is proposed to account for correction terms of increasing order , and its practical implementation and limitations are discussed .
we review previous studies employing this novel approach on an array of ( bio)molecular systems providing a solid basis for its application and demonstrating the importance of these correction terms in the evaluation of absolute entropy and entropy differences . using microsecond md simulation of a test system
, we analyze the convergence properties of the absolute single - molecule entropy and of the correction terms used to estimate it .
the results emphasize that sufficient ps sampling is required for a reliable estimation of entropic contributions because convergence of both the qh upper bound and the required correction terms should be achieved .
qh analysis aims to account for motions in the overall extent of ps accessible to a molecular system at thermodynamic equilibrium .
it relies on approximating the configurational probability distribution as a multivariate gaussian , the momenta of which can be estimated , e.g. , from molecular dynamics ( md ) or monte carlo simulations .
more precisely , for a given choice of generalized coordinate system q ( of dimension m = 3n , n being the number of atoms ) , its input quantity is the covariance matrix cq characterizing the atom - positional fluctuations ( and their correlations ) around an average configuration q. assuming a canonical ensemble and fluctuations resulting from an underlying harmonic potential of the form 

 where 
 is an effective hessian matrix and qo an effective equilibrium configuration , it follows that(11 ) 
 note that the corresponding harmonic model only strictly produces the correct average configuration q and covariance matrix cq for generalized coordinate systems where the mass - metric tensor aq is configuration independent.(11 ) in this study , we only consider the specific case of single - molecule entropy ( i.e. , the entropy of individual distinguishable atoms in a covalently bound molecule ) based on md simulation trajectories . as detailed elsewhere,(21)single - molecule entropy differs from molecular entropy in that the former estimate only accounts for intermolecular correlation in terms of the effect of the solvent on the single - molecule dynamics . in practice ,
the qh analysis of an md trajectory involves the following steps.(11 ) first , the average configuration q and the covariance matrix cq in the chosen coordinate system are evaluated as 

 the equilibrium configuration qo and hessian matrix 
 of the effective underlying harmonic model are then defined according to eq 2 .
second , the ( symmetric ) metric - tensor - weighted covariance matrix is diagonalized 

 where vq is a m  m - dimensional ( orthogonal ) matrix the columns of which represent the m components of the eigenvectors { q , m|  m = 1 , ... , m } ( called qh modes ) of the metric - tensor - weighted covariance matrix and fq is a diagonal matrix containing the corresponding eigenvalues .
these eigenvalues are related to the associated angular frequencies of the underlying effective harmonic model as ( see eqs 2 and 4 ) 

 the sum of the eigenvalues in fq is equal to the total mean - square metric - tensor - weighted fluctuation of the system , i.e. 

 so that the eigenvalues can be interpreted as contributions of individual qh modes to this quantity ( a larger value indicating a larger contribution to the total fluctuation of the molecule ) .
third , the simulated trajectory is projected onto the qh modes , i.e. , one considers the transformed coordinates bq defined as 

 these so - called qh coordinates satisfy the properties(11 ) 
 because fq is diagonal , eq 8 enforces that the individual components { bq , m|  m = 1 , ... , m } of the qh
coordinates bq are pairwise linearly uncorrelated , which , however , does not imply the absence of higher order ( i.e. , pairwise supralinear and higher order ) correlations .
we previously motivated the choice of a cartesian vs internal coordinate system.(11 ) if a cartesian coordinate system r is employed ( after removal of the overall translational and rotational motion from the sampled trajectory(22 ) ) , the mass - metric tensor ar is identical to the mass matrix m ( thus configuration independent , so that eq 2 is exactly satisfied ) . in this case , the qh analysis relies on the diagonalization of the mass - weighted cartesian covariance matrix , i.e. 

 in place of aq  cq  aq in eq 4 . in the absence of geometric constraints ,
if nc geometrical constraints are present in the system ( e.g. , bond - length constraints ) , these will map to an identical number of zero eigenvalues ( see appendix a in ref ( 23 ) for a derivation in the mathematically similar context of essential - mode analysis ) .
thus , the number of qh modes with nonzero eigenvalues is m = 3n  nc  6 , where m = 3n . when using a generalized coordinate system excluding overall translation and rotation variables , one has m = m = 3n  nc  6 .
note that the qh coordinates have units of mass  length.(11 ) single - molecule entropies can be obtained as follows.(11 ) in terms of qh coordinates , the configurational probability distribution associated with the effective harmonic model of eq 2 corresponds to that of m independent harmonic oscillators .
assuming a canonical ensemble and a configuration - independent mass - weighted metric tensor , this leads to(11 ) 
 where s( ) is the canonical entropy of a one - dimensional harmonic oscillator with angular frequency . the classical expression scl , o( ) and the quantum - mechanical expression sqm , o( ) for this quantity are 

 and 

 where  = h(2 ) is the reduced planck s constant , leading to eq 10 to corresponding total estimates scl , o and sqm , o , respectively . in practice , even if the underlying trajectory was generated at the classical level , the qh entropy must be evaluated using the quantum - mechanical oscillator formula because in the high - frequency limit the classical entropy of a one - dimensional harmonic oscillator diverges to the unphysical limit of  rather than to the physical limit of zero .
however , the qh entropy estimate sqm , o is not the absolute configurational entropy of a single molecule but an upper bound for this quantity due to the presence of qh mode anharmonicities and correlations not accounted for in the effective harmonic model of eq 2 .
corresponding correction terms can be formulated exactly at the classical level using an approach previously described(11 ) and briefly summarized below . in the canonical ensemble , assuming a configuration - independent mass - metric tensor , the exact classical single - molecule entropy reads(11 ) 
 where p(bq ) is the probability distribution in the m - dimensional space of the qh coordinates bq ( eq 7 ) .
this expression can be compared with the approximate ( classical ) qh estimate scl , o based on eqs 10 and 11 , i.e. 

 a series of increasingly accurate estimates { scl , k|  k= 0 , 1 , ... , m }
may now be formulated as 

 where c denotes a combination of k qh modes , c(k , m ) = [ ( m  k)!k!]m ! for k > 0 along with c(0,m ) = 0 represents the total number of possible combinations c of k modes among the m qh modes and p(bq ) is the k - dimensional probability distribution in the subspace of the qh coordinates bq within bq that are involved in a combination c. the derivation of this equation is given in the appendix in the supporting information .
it is easily verified that scl , k=0 = scl , o ( eq 14 , i.e. , the uncorrected classical qh entropy ) and scl , k = m = scl ( eq 13 , i.e. , the exact classical entropy ) .
substituting the classical estimate scl , o by the corresponding quantum - mechanical estimate sqm , o ( eqs 10 with 12 ) into eq 15 and introducing successive correction terms defined as 

 leads to a ( classically ) corrected qh entropy estimate 

 the successive correction terms of eq 17 involve integrals over the probability distributions p(bq ) in eq 15 with increasing dimensionality k. note that these terms are all individually negative ( or vanishing ) .
the first correction term scl,1 involves one - dimensional ( 1d ) integrals and accounts for anharmonicities in the individual qh modes .
the second correction term scl,2 involves two - dimensional ( 2d ) integrals and accounts for pairwise ( supralinear ) correlations between the qh modes . for simplicity
, these two terms will be renamed sclah and sclpc , respectively , to match the notation used in other studies .
the following higher order correction terms account for correlations among qh modes beyond the pairwise ones .
although the classical qh entropy estimate scl , o usually represents a poor approximation to its quantum - mechanical counterpart sqm , o , the evaluation of the correction terms at the classical level remains accurate because anharmonicities and correlations principally affect the low - frequency qh modes for which the classical approximation holds.(11 ) the successive correction terms in the series of eq 17 are increasingly difficult to evaluate because both ( i ) the number of terms c(k , m ) involved in the evaluation of scl , k and ( ii ) the sparseness in the required multiple - mode probability distributions p(bq ) increase exponentially with k. for this reason , their evaluation is restricted in practice to the first two terms and implies an intrinsic uncertainty on the final estimate compared to the true single - molecule entropy ( i.e. , persisting in the limit of infinite sampling ) . however , note that , in a different context , alternative approximate formulations to estimate terms of increasing order mutual information have been proposed and seem to suggest that the first two correction terms in eq 17 are indeed dominant .
following from eqs 16 and 17 , the expressions for the first two correction terms are 

 and 

 leading to the corrected absolute single - molecule entropy estimate 

 the relative magnitudes fclcrc , fclah , and fclpc of the correction terms sclcrc , sclah , and sclpc with respect to the qh entropy upper - bound sqm , o ( expressed in percent ) , i.e. 

 may then serve as a measure for the importance of the aforementioned corrections . in practice ,
the 1d and 2d integrals involved in eqs 18 and 19 are evaluated numerically in the form of sums over corresponding histograms .
it is reasonable to choose the bin width along a given qh mode in proportion to the width ( first moment ) of the probability distribution along this mode with proportionality factors 1 and 2 for 1d and 2d integrals , respectively .
however , 1 and 2 values must be selected carefully in order to keep both finite - sampling and binning errors to a minimum , i.e. , to ensure the independence of the results on these two parameters.(11 ) for this reason , we monitored the dependence of such numerical integrals on the width of histogram bins for increasing periods of time , as described in section .
note , finally , that the absolute single - molecule entropies so far discussed exclude roto - translational contributions . in principle
, a translational entropy contribution can be included using the quantum - mechanical expression of the sackurtetrode equation for a specified standard state of the pressure ( molecule in the gas phase ) or of the concentration ( molecule in solution ) .
similarly , the rotational entropy contribution could be included using the appropriate quantum - mechanical expression ( e.g. , rigid - rotor approximation , based on the average inertia tensor of the molecule ) . however , these two contributions are likely to be highly coupled with each other and with s , i.e. , they are not strictly additive , and their rigorous treatment is therefore still challenging . a recent study reported on relatively small effects of motional correlation on changes of reorientational entropy using selected qh modes from a 1.5 ns simulation of the ubiquitin protein.(35 ) in the present article , single - molecule configurational entropies refer to entropies excluding roto - translational effects .
a 1.1 s long md simulation of the cc peptide ( ch3-co - s - i - r - e - l - e - a - r - i - r - e - l - e - l - r - i - coo ) at 300 k was performed with the amber 9 software,(36 ) the amber 99sb parameter set,(37 ) and the compatible tip3p water model.(38 ) the simulation was initialized from the -helical configuration based on a x - ray model structure of the cc coiled coil ( pdb i d 1s9z).(39 ) trajectory snapshots were saved every 10 ps for analysis .
the simulation setup and trajectory analyses are detailed elsewhere.(40 ) backbone atom - positional root - mean - square deviations ( rmsd ) from the initial folded structure and radius of gyration ( rgyr ) were calculated using all c atoms .
independent qh analyses were performed for 22 increasingly long segments of the simulation ( differing in length by 50 ns ) by calculation of the solute all - atom mass - weighted covariance matrix dr ( eq 9 ) in cartesian coordinates after least - squares fit superposition(22 ) of all configurations onto the initial structure to eliminate overall translation and rotation and diagonalization ( eq 4 with aq1/2  cq  aq1/2 = dr ) .
a total of 534 ( m = 3  297  351  6 ) modes associated with nonvanishing eigenvalues were considered .
after determination of the qh modes ( columns of the matrix vr in eq 4 ; sorted in order of decreasing eigenvalues , i.e. , increasing m frequency in eq 5 ) , the trajectory was projected in this basis set to obtain the time series of the corresponding qh coordinates br ( eq 7 ) .
this first part of the analysis was performed using the s_correction program as implemented in the gromos++ module of the gromos05 software(41 ) for biomolecular simulation .
the qh entropy upper bound , sqm , o ( eq 10 with eq 12 ) , the corrections for mode anharmonicity , sclah ( eq 18 ) , and pairwise supralinear mode correlation , sclpc ( eq 19 ) , their sum , sclcrc ( eq 20 ) , the improved absolute single - molecule entropy , s ( eq 20 ) , the relative terms fclcrc , fclah , and fclpc ( eq 21 ) , and the sum of the eigenvalues , tr[fr ] ( eq 6 ) , were then calculated for each of the 22 trajectory segments .
note that this analysis is computationally intensive because , as discussed in section , each of the 22 sclah values requires the estimation of 534 1d integrals , while each of the 22 sclpc values requires the estimation of 142  311 2d integrals ( eq 15 ) .
in addition , for each of these integrals the optimized proportionality factors 1 and 2 were determined based on multiple integral calculations for an accurate numerical integration ( see section ) . as an indication of the actual computational cost , using an intel xeon x5450 3.0 ghz , dedicated software , and the above procedure , each
1d or 2d integral can be estimated with an average cpu time of 0.07 or 0.21 s , respectively , from 50 ns trajectory windows .
overall , the analyses presented in this work require a cpu time that sums up to 6 months .
the key findings of previous studies concerning the uncorrected qh upper bound , sqm , o ( eq 10 with eq 12 ) , the improved absolute single - molecule entropy s ( eq 20 ) , and the relative magnitude of the cumulative correction term fclcrc ( eq 21 ) are summarized graphically in figure 1 .
these results span systems with different chemical nature : 2 -peptides in methanol,(11 ) the 11 disaccharides of gluscose in water,(28 ) the dipalmitoylphosphatidylcholine ( dppc ) lipid in a hydrated bilayer,(24 ) the w191 g mutant cavity and its gating loop within cytochrome c peroxidase in water,(20 ) the h - ras lipopetide anchor in water or inserted into a model lipid membrane,(27 ) and the cc-peptide in water ( this study ) . in some cases
, the qh analysis was also performed separately for different chemical environments or conformational states of the molecule , which permits estimating relative entropies , thereby quantifying the impact of the correction terms on the thermodynamic process of interest .
these processes include reversible peptide folding ( ref ( 11 ) and this study ) , conformational changes in carbohydrates(28 ) and lipids,(24 ) lipopeptide insertion in a model membrane bilayer,(27 ) ligand binding to a protein cavity,(20 ) and protein - loop gating.(20 ) the results presented in figure 1 are scaled by the number , n , of atoms to allow for a comparison among molecules of varying size ( the raw data is available as supporting information , table s1 ) .
some clear qualitative trends are evident , although a direct comparison among these studies is not possible due to the different md time scales and physicochemical conditions . in all systems the cumulative correction term sclcrc ( eq 20 ) is generally sizable , demonstrating an overall large deviation from a qh behavior as evaluated up to the pairwise supralinear level . the corresponding relative magnitudes , fclcrc , display values from 9% to 73% of the qh upper - bound value sqm , o ( figure 1 ) .
in detail , these important cumulative terms result from the sum of correction terms for mode anharmonicity ( sclah ; eq 18 ) that are always relatively small ( up to 3% of the upper - bound value sqm , o ) and for pairwise supralinear mode correlation ( sclpc ; eq 19 ) that are always dominant .
the latter correction term has a magnitude that depends on the physical nature of the molecular motional correlation experienced by the molecular system in a given thermodynamic state .
the largest relative corrections , fclcrc , are expected and found for intrinsically more ordered systems ( figure 1 ) .
this can be explained by considering that restricted flexibility is typically promoted by inter- and/or intramolecular interactions , simultaneously inducing increased motional correlation .
for example , the ligand - bound state of the w191 g protein cavity(20 ) displays the largest fclcrc value ( 73% ) , i.e. , the thermodynamic ensemble involving the largest motional correlations and lowest entropy content among those studied . on the other end of the spectrum and in line with this qualitative picture , the smallest fclcrc values were reported for the dppc lipid in a bilayer ( 9%),(24 ) i.e. , the ensemble characterized by the highest molecular flexibility and thus the lowest motional correlations .
interestingly , the 11 disaccharides of glucose in water(28 ) display high variability and always large fclcrc values ( 4572% ) .
this behavior can be explained considering that these molecules involve a reduced number of degrees of freedom overall and the linkage between rather stiff glucose rings is the torsion defining major conformational changes .
these qualitative trends are also in agreement with the observation that entropy is the measure of ps sampling for a molecular system .
the qh upper bound , sqm , o , and the improved absolute single - molecule entropy , s , are estimated based on the ps that has been accessed during a md simulation of finite time scale , i.e. , only a fraction of the ps accessible to the system . these time scales ranged from 50 ns ( w191 g mutant(20 ) showing the lowest entropy ) to 25.6 s ( concatenated trajectory of the dppc lipid,(24 ) showing the largest entropy ) .
however , this hampers a quantitative comparison of sqm , o , s , and sclcrc values among previous studies . prompted by these observations , the dependence of these quantities on the extent of accessed ps was assessed on the microsecond time scale for a peptide test system .
the cc peptide in water was chosen as a test system to investigate entropy convergence properties because of its small size and broad ps accessibility .
figure 2a shows the time series of the backbone atom - positional root - mean - square deviation ( rmsd ) from the folded structure and of the backbone radius of gyration ( rgyr ) along a 1.1 s of md simulation .
the peptide undergoes several reversible folding / unfolding events and samples a variety of unfolded configurations and compact folds.(40 ) cc peptide dynamics on the microseconds time scale and entropy convergence .
( a ) the backbone atom - positional root - mean - square deviation ( rmsd ; black ) from the initial helical fold and of the backbone radius of gyration ( rgyr ; gray ) are shown along the time , t. the cartoon representations highlight example configurations ( oriented with the ch3co terminus down ) .
( b ) build - up curves of the qh entropy upper bound sqm , o ( eqs 10 with 12 ; dashed line ) and of the improved absolute single - molecule entropy s ( eq 20 ; solid line ) .
convergence of ( c ) the cumulative correction term sclcrc ( eq 20 ) and its contribution to the free energy tsclcrc , ( d ) its relative value fclcrc ( eq 21 ) , and ( e ) the sum of the eigenvalues tr[fr ] ( eq 6 ) .
the probability distributions p(br , m ) of the transformed qh coordinates br ( eq 7 ) along selected qh modes ( m = 1 , 2 , 6 , 10 , 50 , and 500 ) are shown in figure 3a .
the actual distributions become increasingly narrow and similar to the gaussian functions for higher m indices , i.e. , the corresponding qh modes become increasingly stiff and harmonic . however , the distributions along the lowest frequency modes ( e.g. , figure 3a , m = 1 or 2 ) differ significantly from gaussian functions and evidently result from the superposition of multiple off - center gaussian - like distributions .
a similar observation was previously reported in the context of two -peptides in methanol for which two main subensembles of folded and unfolded configurations could be disentangled based on the lowest frequency modes ( figures 35 in ref ( 11 ) ) . in the present case ,
the most pronounced peaks for the cc peptide arise from folded configurations ( see figure 3a , m = 1 and 2 , leftmost peak ) .
the actual distributions ( gray line ) are displayed together with the corresponding gaussians , i.e. , po , m(br , m ) = ( 2fr , m)e  br , m the time dependence of these results was investigated as summarized in figure 3b for the two coordinates br with lowest frequencies , i.e. , those contributing the most to the total mean - square metric - tensor - weighted fluctuation of the system ( eq 6 ) .
the corresponding distributions vary significantly with the extent of ps sampling , as revealed by averaging over the first 0.2 , 0.4 , and 0.8 s periods , or over the entire 1.1 s ensemble ( cf .
figure 3b vs figure 3a for m = 1 and 2 ) . increasing the simulation time results into broader distributions due to the larger extent of ps sampled .
the intensity of the leftmost peak , corresponding to the contribution of the folded configurations , clearly reduces along the simulation initialized from the cc helical fold and evolving through a broad range of heterogeneous configurations ( figure 1a ) .
the data indicate that convergence of the probability distributions associated with the low - frequency qh coordinates in br requires sampling times longer than 1 s , considering that such differences persist when comparing results from the first 800 ns period with the whole 1.1 s simulation . the entropy convergence for the cc peptide in water as a function of sampling time is illustrated in figure 2b .
the upper - bound curve obtained by application of the uncorrected qh formula ( sqm , o ; eq 10 with eq 12 ) is compared to the build - up curve of the improved absolute single - molecule entropy ( s ; eq 20 ) .
interestingly , the qh upper bound reaches convergence noticeably faster than the improved absolute single - molecule entropy . in detail , 0.3 or 0.7 s
is needed to sample 90% or 99% of the final sqm , o estimate of 6922 j k mol , while larger sampling times of 0.5 or 1.0 s are needed to sample 90% or 99% of the final s estimate of 5916 j k mol ( see also supporting information , table s2 ) .
these results clearly demonstrate that the convergence of the sqm , o upper bound does not imply the convergence of the absolute single - molecule entropy , s. in fact , the first quantity relies on the convergence of linear motional correlations only ( following from the definition of linearly independent qh modes in eq 7 with eq 8) . instead ,
as described below , the second quantity requires in addition the convergence of supralinear motional correlations . for this reason
, s seems to represent a better indicator of convergence ( compared to sqm , o ) for the absolute single - molecule entropy .
the convergence behavior of the cumulative correction term , sclcrc ( eq 20 ) , was also monitored ( figure 2c ) .
this entropy term describes the overall deviation from the qh model due to mode anharmonicity ( sclah , eq 18 ; 4 j k mol after 1.1 s ) and correlation ( sclpc , eq 19 ; 1002 j k mol after 1.1 s ) effects , associated with all unique combinations of modes m and n. its convergence behavior can be used as well as a measure of uncertainty on the entropy estimate . for comparison with previous studies ( figure 1 ) , the time evolution of the corresponding relative contribution fclcrc ( eq 21 ) to the uncorrected sqm , o estimate ( eq 10 with eq 12 ) is also displayed ( figure 2d ) . in all studies , the dominant part of this correction arises from qh pairwise ( supralinear ) mode correlations ( fclah values are < 0.05% ; see also supporting information , tables s1 and s2 ) .
importantly , it is found that the magnitude of the correction term sclcrc monotonically decreases from an initial value of 2745 j k mol ( first 0.5 s sampling ) to a final value of 1006 j k mol ( entire 1.1 s sampling ) , showing an initial convergence behavior .
this suggests that limited ps sampling results in both the underestimation of sqm , o and the overestimation of sclcrc ( predominantly through its sclpc component ) , both artifacts leading to an underestimation of the final absolute single - molecule entropy s. this result can be explained considering that motional correlations are larger for a molecular system sampling a confined part of ps as opposed to sampling of a multiple - minima landscape . for the cc test system
these results demonstrate that a limited ps sampling leads to the overestimation of corresponding motional correlations , thus of sclpc values with respect to that expected for a canonical ensemble of the same system at thermodynamic equilibrium .
the mass - weigthed root - mean - squared fluctuation , i.e. , the sum of the eigenvalues of the mass - weighted covariance matrix , tr[fr ] ( eq 9 ) , was also monitored along time as an independent measure of convergence ( figure 2e ) . in terms of tr[fr ] , we note that a first plateau region is reached after 1.1 s ( figure 2e ) , in line with the s values ( figure 2b ) .
this observation confirms as well that the cc peptide is not trapped in a few local minima .
instead , it explores new configurations even after several hundreds of nanoseconds ( figure 2a ) .
this is evident when the term is expressed in the form of its contribution to the system free energy , tsclcrc ( figure 2c , right axis ) .
the resulting value ( 302 kj mol based on 1.1 s ) is about an order of magnitude larger than the free energy changes of typical ( bio)chemical processes .
thus , although partial cancellation of this term can be expected for entropy differences between two different molecular environments or conformational states , small differences will still lead to large free - energy contributions ( of sign and magnitude difficult to be predicted a priori ) .
in addition , this result suggests that time convergence of the entropy estimate should be taken into account as well when comparing the efficiency and accuracy of alternative computational approaches .
second , we stress that all m per - mode contributions need to be included for an accurate estimation of eqs 10 and 18 because modes with large m indices ( high frequencies ) also contribute to sqm , o ( data not shown ; see figure 8 in ref ( 11 ) for a similar analysis ) .
this marks a difference with what is typically observed for the contribution of a reduced number of essential modes to the total system fluctuation.(23 ) due to the similar mathematical formalism,(11 ) this argument can be easily demonstrated as well for the calculation of entropies from normal - mode analysis for systems sampling one local ps minimum .
third , the analysis of the leading correction sclpc in terms of all c(2,m ) = [ ( m  2)!2!]m ! unique pair combinations reveals that not only modes with low indexes ( high amplitudes , low frequencies ) contribute substantially .
thus , all pairs of qh modes need to be considered in eqs 15 and 19 .
this requirement arises from the observation that high correlations can be present among modes with either large or small m and n indices ( low or high frequencies ; figure 12 in ref ( 11 ) ) .
interestingly , this behavior was observed for highly flexible systems ( e.g. , the cc peptide of this study or the reversibly folding -peptides in ref ( 11 ) ) but not for more rigid systems confined to a local ps sampling ( e.g. , the w191 g cavity in ref ( 20 ) , unpublished results ) . whether the latter result depends on a limited ps accessed or on the physical nature of qh - mode correlation remains to be addressed .
the analysis presented in this work relies on the numerical integration of the actual probability distributions p(br , m ) and p(br , m , br , n ) evaluated based on the md trajectory ( eqs 18 and 19 ) .
two alternative procedures were described to estimate these 1d and 2d integrals with optimal ( nonarbitrary ) histogram bin widths , as detailed in appendix c of ref ( 11 ) . in this study ,
optimal parameters 1 and 2 were chosen as the midpoint between the intersections of a horizontal line with the limiting lines for too fine and too coarse integration at the optimal value of the 1d or 2d integrals in the graph showing these values as a function of ln 1 or ln 2 , as summarized in figure 4 .
dependence of the numerical integration of probability distributions on the width of histogram bins for increasing periods of time , t. ( a ) integrals over the 1d distributions involved in eq 18 are shown for eigenvectors 1 ( circles ) , 50 ( squares ) , and 500 ( crosses ) .
( b ) integrals over the 2d distributions involved in eq 19 are shown for eigenvector pairs 1,2 ( circles ) , 1,100 ( squares ) , and 1,500 ( crosses ) .
the results are displayed for the whole ensemble ( 1.1 s ; black ) or t = 200 ( red ) , 400 ( green ) , 600 ( blue ) , 800 ( cyan ) ns as a function of ln 1 ( a ) or ln 2 ( b ) , where  is the ratio of the bin width along each dimension to the corresponding distribution width .
the middle point between a pair of limiting lines for too fine ( left side ) and too coarse ( right side ) numerical integrations ( dashed lines ) defines optimal 1 or 2 values .
figure 4a shows the values of the 1d integrals for a sample set of eigenvectors ( m = 1 , 50 , and 500 ) , evaluated numerically using different values of 1 .
values approaching these limiting curves are incorrect because they show a dependence of the evaluated integral on the bin size .
however , for each curve , a clear plateau defines the range of 1 values for which the integration result is essentially independent of the bin size .
finite - sampling artifacts affect the integration with the smallest values of 1 , while coarse - binning artifacts affect the integration with the largest values .
figure 4b shows the values of the 2d integrals for a sample group of eigenvector pairs ( m , n = 1,2 ; 1,100 ; 1,500 ) , evaluated numerically with different values of 2 , together with the corresponding limiting lines and the optimal 2 values . here
, the plateau regions are narrower and the value of 2 has to be chosen more carefully .
note that 2d integrals are always negative when estimated using the optimal 2 values , but incorrect positive values would be obtained based on too small 2 values .
the dependence of these curves on the simulation time was also monitored ( figure 4 ) .
all 1d and 2d curves show a coarse - integration limit that is essentially independent of the md period considered . yet
, they also show that the fine - integration - limiting curves shift to lower 1 and 2 values upon increasing the simulation time , thus reducing the dependence on the integration bin size .
this effect is more pronounced for the 2d integrals because they require more data points than 1d integrals .
overall , these results demonstrate that the procedure employed in this work allows estimating both 1d and 2d integrals of eqs 18 and 19 in a nonarbitrary way as a function of the simulation time .
the presence of plateau regions independent of the integration bin size for all md periods considered shows that the observed change of integral values along the simulation time largely depends on the extent of ps sampled yet not on the numerical procedure employed . a similar analysis performed in the case of corresponding 3d integrals ( triplewise combinations of qh modes )
revealed that indeed no such behavior can be achieved , although using a 1.1 s trajectory . in practice ,
eq 15 can be only estimated for the first two terms owing to finite sampling artifacts and data sparseness .
the theory and practical implementation of an approach recently proposed(11 ) to estimate improved configurational entropies from quasi - harmonic analysis of molecular dynamics simulations are briefly reviewed .
it involves the calculation of correction terms of increasingly high order to account for deviations from the quasi - harmonic approximation in frustrated molecular systems .
the convergence properties of the absolute single - molecule entropy are critically investigated using microsecond molecular dynamics simulation of the cc peptide in water .
prompted by the comparison of the results with previous studies addressing mode anharmonicity and correlation effects , the convergence behavior of individual quasi - harmonic modes , of the absolute single - molecule entropy , and of the correction terms for anharmonicity and pairwise ( supralinear ) correlations are analyzed .
our data provide a number of new insights to tackle the challenge of accurate entropy estimation by computer simulation . in line with a previous study,(11 )
the probability distributions associated with components of the quasi - harmonic coordinates only deviate significantly from gaussian functions for the first few components , resembling the behavior observed in the different context of a single - atom displacement for an -helical peptide(42 ) and of essential modes for protein dynamics.(23 ) for these components , the probability distributions result from a superposition of clearly distinguishable contributions from the folded and unfolded ensembles .
however , it is shown that the components of these eigenvectors converge slowly ( > 1 s ) , consistent with the observation that the cc peptide steadily explores new configurations . in line with previous studies ,
the entropic contribution of anharmonicity is small while the pairwise ( supralinear ) correlation correction to the entropy is large .
the deviation from the quasi - harmonic assumption affects more significantly conformational states dominated by high motional correlation . using microsecond molecular dynamics simulation of a peptide test system
we show that limited phase - space sampling results in an overestimation of correlation effects , and we discuss its implications for entropy estimation .
this study demonstrates that the convergence of the quasi - harmonic upper - bound entropy with simulation time does not imply the convergence of the system absolute single - molecule entropy . as a consequence
, our study also suggests that the convergence of the absolute single - molecule entropy rather than that of the quasi - harmonic upper bound should be preferably monitored .
because the cumulative correction term accounting for both mode anharmonicity and pairwise ( supralinear ) correlation effects converges slowly and monotonically decreases , previous studies based on shorter time scales may have , in some cases , partly overestimated this correction term , thus leading to underestimated absolute entropy estimates .
overall , the present study emphasizes the need of sufficient phase - space sampling to estimate entropic contributions from computer simulations .
ideally , only thermodynamic ensembles at equilibrium should be considered to this end , i.e. , full phase - space sampling obtained from simulations on time scales of several microseconds . in practice
, we suggest that enhanced sampling techniques and/or concatenated copies of independent simulation trajectories will be useful tools to alleviate these problems in the future if properly combined with the correction terms used herein.(28 ) this strategy will open the possibility to include as well correlation effects of higher order than the pairwise ( supralinear ) explicitly considered in this study .
a bright future opens for the estimation of accurate thermodynamic properties for biomolecular systems using chemical theory and computation .
 the convergence properties of the absolute single - molecule configurational entropy and the correction terms used to estimate it are investigated using microsecond molecular dynamics simulation of a peptide test system and an improved methodology . 
 the results are compared with previous applications for systems of diverse chemical nature . 
 it is shown that ( i ) the effect of anharmonicity is small , ( ii ) the effect of pairwise correlation is typically large , and ( iii ) the latter affects to a larger extent the entropy estimate of thermodynamic states characterized by a higher motional correlation . 
 the causes of such deviations from a quasi - harmonic behavior are explained . 
 this improved approach provides entropies also for molecular systems undergoing conformational transitions and characterized by highly frustrated energy surfaces , thus not limited to systems sampling a single quasi - harmonic basin . 
 overall , this study emphasizes the need for extensive phase - space sampling in order to obtain a reliable estimation of entropic contributions . 
Introduction
Methods
Results and Discussion
Conclusion