equilibrium constants 
 ( and their derivatives ) that describe the 
 delicate balance between the native and denatured states of proteins 
 are frequently measured to understand if and why certain proteins 
 are more stable than others and how mutations can alter this stability . however , under physiological conditions , the 
 ratio of [ denatured]/[native ] is usually very small , making it difficult 
 to measure accurately . often the equilibrium must be perturbed , for 
 example by changing the temperature or pressure , in order to achieve 
 an equilibrium constant that is larger in magnitude and can , therefore , 
 be calculated more precisely . afterward , the results are usually extrapolated 
 back to physiological or standard conditions . thus , it is important to understand the effects of temperature and 
 pressure on the equilibrium constant . however , structural properties are usually the sole or primary focus 
 of most simulations studies , and the corresponding thermodynamic properties 
 are often ignored . when thermodynamic properties 
 are determined from simulation , it is typically by way of a method 
 that does not retain information concerning which conformation gave 
 rise to the specific values , meaning the simulation does not really 
 provide any more insight than an experiment would , or the properties 
 are calculated from phenomenological relationships rather than rigorously 
 derived thermodynamic expressions . although multiple environmental variables affect 
 the conformational 
 equilibrium of a protein , here we limit 
 our scope to thermal and pressure effects . hence , the effects of cosolvents 
 and ph on protein equilibrium are not discussed here , although a theoretical 
 framework for understanding their effects on an equilibrium has recently 
 been advanced . to understand the effects of the environment on a protein s 
 thermodynamic properties , it is common to study an infinitely dilute 
 protein in a single solvent . to keep the situation tractable , we will 
 assume that a protein molecule adopts one of only two forms , the native 
 state , n , or the denatured state , d. ( the process of classifying a 
 specific conformation as being n or d is subjective , but addressing 
 this subjectivity is beyond the focus of this paper . when studying 
 a specific system s equilibrium by computer simulation , it 
 would be wise to use the same definition as was used experimentally , 
 e.g. , measuring tryptophan fluorescence , if possible . ) when considering 
 the equilibrium between n and d , n  d , the standard state 
 gibbs free energy of unfolding ( ndg , hereafter g ) determines the population of the protein conformations 
 through the usual relationship to the thermodynamic equilibrium constant k = ad / an , where ai is 
 the equilibrium activity of form i , through1where r is 
 the gas constant , t is the temperature , i is the chemical potential of the protein 
 in conformation i , and the superscripted o denotes 
 the standard state . the activity coefficients 
 at infinite dilution are assumed to be one , which leads to the replacement 
 of the activities with their concentrations in the equilibrium constant 
 expression . while g quantifies the ratio of species 
 at equilibrium , derivatives of g reveal 
 increasingly more detailed information the higher their order . the 
 bivariate taylor expansion of the standard state free energy change 
 for unfolding , g( , p ) = ln(k / kref ) , up to second order around a reference ( ref ) inverse temperature 
 (  = 1/rt ) and pressure ( p ) provides an avenue to interpret experimental thermal and pressure 
 denaturation data via2where k is the column vector 
 whose elements are the first derivatives of ln k with 
 respect to  and/or p evaluated at a reference 
  and p and k is the matrix whose 
 elements are the corresponding second partial derivatives of ln k. d is the column vector of deviations from 
 the reference  and p,3where x = x  xref . it can be shown 
 that4where h denotes 
 the enthalpy , v the volume , cp the heat capacity , p the thermal expansion , and t the isothermal compressibility . the changes are given by h  hd  hn , v  vd  vn , cp  cp , d  cp , n , p  vdp , d vnt , n. first order taylor 
 expansions provide the following temperature 
 and pressure dependencies of h , the 
 standard state entropy change ( s ) , and 
 v via,5thus , to use md simulations 
 to understand how temperature and pressure affect the conformational 
 equilibrium of a protein in significant detail , that is , to map g as a function of  and p , 
 one requires h , v , cp , p , and t. similarly , to map h , s , and v as a function of  and p , 
 one requires cp , p , and t. the second derivative properties 
 can also be expressed by6the relationships in eq 6 will be used later . 
 all of the above properties are standard state 
 properties . standard state properties are usually associated with 
 those of an infinitely dilute protein ( for which the activity is approximated 
 by the concentration ) in either the native or denatured forms , which is helpful since most molecular dynamics 
 simulations model the infinitely dilute case . in fact , the single 
 protein molecules studied by simulation should really be considered 
 as pseudo infinitely dilute ( but we will drop the 
 use of pseudo for the sake of brevity ) , since the amount of protein 
 in a given volume ( the formal or effective concentration ) will vary 
 depending upon the total system size even though the number of protein 
 molecules is fixed at one . most simulations do not readily provide 
 the above infinitely dilute 
 partial molar properties , because the expressions relating the output 
 simulation results to these properties have generally been considered 
 unknown . in contrast , it is well - known that the system fluctuations 
 are related to the compressibility , thermal expansion , and heat capacity 
 of the whole system . in the npt ensemble the corresponding 
 expressions are7where the angled brackets 
 indicate time or 
 ensemble averages , x = x  
 x and xy = xy we emphasize that these are all bulk system properties , 
 which means that when applied to a system containing an infinitely 
 dilute solute , the properties obtained are those of the solvent , perturbed 
 to some extent by the presence of the solute molecule . for a binary 
 mixture in which both species 1 and 2 are at finite concentrations , 
 the process of extracting partial molar properties is straightforward . one needs only to fit the property of interest as a function of composition 
 and then take the partial derivatives with respect to each component . however , for the infinitely dilute case , an 
 appropriate method for determining partial molar properties has not 
 been clear or computationally intuitive . computationally speaking , 
 we are only aware of one method in the literature that has been used 
 to obtain all of the infinitely dilute partial molar properties and 
 can be considered thermodynamically rigorous , namely , the enhanced 
 sampling technique of replica exchange molecular dynamics ( remd ) . remd simulations map out the phase space and 
 allow one to obtain k by determining the fraction 
 of protein molecules that are folded or unfolded as a function of 
 temperature . once this is known , the free energy of unfolding as a 
 function of temperature can be computed according to8where i is the number density of i ( ni/v ) . after 
 the free energy of unfolding has been computed , the difference ( unfolded 
 minus folded ) in all of the thermodynamic properties may be fitted 
 simultaneously according to the approach illustrated by , for example , 
 garcia et al . note , however , that all 
 that is obtained using remd is the difference between the thermodynamic 
 properties in the n and d forms , not the values of the thermodynamic 
 properties for each of the native and the many denatured states that 
 collectively gave rise to that difference . thus , the remd approach 
 provides an objective way to calculate the change in an infinitely 
 dilute partial molar quantity , x2 , and is 
 therefore extremely useful from a force field validation or a property 
 prediction point of view . what remains missing , however , is the ability 
 to assign thermodynamic properties to specific conformations , that 
 is , the ability to rank different conformations according to their x2 values . many approaches for the determination of partial 
 molar quantities of individual biomolecule conformations have been 
 proposed in the literature . from a computational 
 viewpoint , many subjective definitions for the volume of a protein 
 exist . levy and co - workers have used three possible definitions to 
 calculate the volume of 15 proteins . the 
 most striking volume variability was , naturally , for the smallest 
 protein they studied , the insulin monomer ( 51 residues ) . they reported 
 the van der waals ( vdw ) volume to be 5.564 nm , the molecular 
 volume ( calculated using a sphere with probe radius of 0.14 nm ) to 
 be 6.848 nm , and the excluded volume ( also calculated 
 using a sphere with probe radius of 0.14 nm ) to be 11.142 nm . this corresponds to a 100% increase 
 in the volume on going from a vdw volume definition to an excluded 
 volume definition . even the largest protein they studied , carboxypeptidase 
 a ( 307 residues ) , had a 17% increase in the volume on going from a 
 vdw to an excluded volume definition ( reported volumes : vdw = 33.571 
 nm , molecular = 42.098 nm , excluded = 58.055 
 nm ) . while there is nothing 
 wrong with calculating these quantities and using them to probe the 
 properties of proteins , it is unclear which approach best corresponds 
 to the rigorously defined partial molar volume . unfortunately , for 
 those interested in further calculations that depend upon the protein 
 volume , it is also unclear which volume definition should then be 
 used . we have found several examples where an objective definition of the protein volume has been used 
 via the apparent molar approach usually encountered in experimental 
 studies , that is,9where v is 
 the average system volume from an npt simulation , v1 * is the molar volume of pure 1 , n1 is 
 the number of molecules of species 1 , and the angled brackets , again , 
 refer to time or ensemble averages . buff 
 integral , g21 , is given by11and where g21(r ) is the radial distribution function 
 of the solvent ( 1 ) 
 around the solute ( 2 ) defined in the grand canonical ( vt ) ensemble . equation 10 has been 
 used by several authors . due to the nonspherical geometry of proteins , the 
 integration in eq 11 is more difficult to perform 
 if one is interested in a surface based approach . similarly , hirata and colleagues have derived statistical 
 mechanical expressions for the infinitely dilute partial molar volumes 
 using the reference interaction site model ( rism ) integral equation 
 theory coupled to kb theory . it is worth noting 
 that while traditional kb theory provides a 
 rigorous expression for the volume only , the generalization of kb 
 theory to fluctuation solution theory ( fst ) provides rigorous expressions 
 for all of the thermodynamic properties considered in this work . however , 
 fst is still in its infancy , and therefore , while the expressions 
 have appeared , they have not yet been applied to obtain the full set 
 of possible thermodynamic properties . most 
 calculations of t,2 we have found 
 can not be considered rigorous for one of two main reasons . first , 
 they may have relied on subjective definitions of the protein volume . second , they often assume that the equation that relates the bulk 
 system compressibility to fluctuations in the volume of the entire 
 system can be applied to a component of the system ( i.e. , the protein ) 
 to obtain its compressibility , that is , that t,2 = (v2)/v2 , where v2 is some measure 
 of the protein volume . this may or may not provide a measure of the 
 intrinsic compressibility of the protein , but it is unlikely that 
 it corresponds to the thermodynamically defined infinitely dilute 
 partial molar compressibility of the protein . this idea has been used 
 in the literature many times and appears to have originated from the 
 work of cooper . it has been further supported 
 by hinz and co - workers , but has also been strongly 
 contested by imai and hirata . furthermore , 
 while many researchers have been able to successfully obtain reasonable 
 agreement with the experimental t,2 by 
 breaking the simulated compressibility into different contributions , 
 such as the intrinsic ,   hydrational ,  
 and/or  thermal  components ( which involves the use 
 of nonunique volume definitions ) , here 
 we seek to calculate the thermodynamically defined property , t,2. a few noteworthy exceptions to the above compressibility 
 studies 
 exist , such as the work of yu and co - workers mentioned in the previous 
 paragraphs explaining protein volume calculation methods . yu used 
 eq 10 and then looked at the pressure dependence 
 of the volume to obtain the corresponding compressibility . likewise , hirata and co - workers also used the 
 pressure derivative of the volume that was obtained from rism combined 
 with kb theory to obtain the isothermal compressibility . in contrast to these studies , we seek a method that can be used 
 to analyze computer simulations , does not involve integral equations , 
 and directly provides any first or second derivative of the free energy , 
 rather than only providing the volume ( obtained from the first pressure 
 derivative of the free energy ) . in a computer simulation , the absolute partial molar enthalpy is 
 available , whereas experiments only provide the excess partial molar 
 enthalpy . when a single protein conformation is simulated , it is not 
 obvious that there would be any value in quoting the absolute enthalpy 
 of that protein form this may explain why we have found considerably 
 fewer attempts to extract h2 or cp,2 from simulations as compared to the preceding properties , 
 outside of the remd studies . for interesting explanations as to the 
 dearth of simulations targeting cp,2 , the reader is referred to a review by prabhu and sharp and an article by cooper . jones ( lj ) sphere in a lj binary mixture by 
 the finite difference of the system internal energy in the absence 
 or presence of the infinitely dilute solute . however , the method was not generalized to all properties or other 
 types of systems such as biomolecules . as we will show here , it is 
 trivial to make these generalizations , but the literature from the 
 intervening years suggests that this approach has not been popular . in the same vein , smit and co - workers recently calculated the enthalpy 
 of a lipid bilayer , with respect to the enthalpy of the solvent background , 
 in order to make a connection to differential scanning calorimetry 
 experiments . . therefore , smit calculated 
 the ( extensive ) enthalpy of a lipid aggregate , not an ( intensive ) 
 infinitely dilute partial molar lipid enthalpy . however , their method 
 is in the same spirit as the method used in this work and it could 
 be applied to any infinitely dilute solute ( considering the lipid 
 aggregate as a single object ) and to any property , not just the enthalpy . lastly , very few examples of protein thermal expansion values calculated 
 from computer simulations studies can be found . they have all used 
 nonunique protein volume definitions and/or decomposed the thermal 
 expansion into different contributions , such as the intrinsic ,  
  hydrational ,  and/or  thermal  components , 
 which introduces several additional subjective parameters . the purpose of 
 this paper is to investigate 
 the computational feasibility of using a general method for extracting 
 the thermodynamically defined x2 that does 
 not result in a loss of information about the conformational state 
 of the protein that gave rise to these properties . the significance 
 of this work is that it allows md simulations to better fulfill their 
 claim of providing exquisite details about systems , often unavailable 
 from experiment , and to do so using thermodynamically exact relationships 
 and approaches . as mentioned in the introduction , obtaining thermodynamically rigorous partial molar properties for 
 species whose concentrations are finite is trivial by computer simulation . 
 however , finite concentrations are typically not feasible for biological 
 system simulations . thus , it is the infinitely dilute limit that is 
 of primary interest for this work . the current approaches found in 
 the literature to determine the relevant properties are often reasonable , 
 but are typically subjective and often rather involved . in comparison , 
 we show here that it is often simpler to use a thermodynamically rigorous 
 apparent molar approach . as discussed in the introduction , obtaining 
 partial molar properties is , in principle , straightforward when all 
 species are at finite concentrations . for a multicomponent system 
 of 1  solvent , 2  solute , 3 , ... , n  cosolvents , a molar property is defined as xm = x / n , where x is the corresponding extensive thermodynamic property 
 of interest and n is the total number of molecules 
 in the system . here we will be focused on x = h , v , cp , vt , and vp only , where we have multiplied 
 the t and p by the system volume . in addition , we have also investigated 
 results for the kinetic energy , k , but only as a 
 test of the procedure , so as to ensure that the equipartition results 
 were achieved . the partial molar property of 2 is then obtained from12where the subscripted { n} 
 indicates that all n other than n2 are held constant . generally , 
 it would not be feasible to perform simulations in which one varies 
 the number of protein molecules , because the required system sizes 
 would become elephantine rather quickly . furthermore , researchers 
 are often more interested in the infinitely dilute system where the 
 inclusion of protein  protein interactions would often be undesirable . 
 a finite difference approach can be used if n2 is set equal to zero and to one in two separate simulations 
 while { n} is held fixed . as mentioned in 
 the introduction , this finite difference calculation 
 has been used before for calculations of a lj sphere s volume 
 and internal energy , for determining protein volumes , and for a lipid 
 aggregate enthalpies . it has also been noted that such an approach can suffer from large 
 statistical errors . we instead 
 have constructed a series of infinitely dilute systems 
 in which we vary x2 by changing n1 , where n2 is always 
 one . if we choose systems such that even the smallest system contains 
 water that is beyond the  sphere of influence  of the 
 protein , that is , some waters exhibit the thermodynamic properties 
 of bulk water , then the addition of more waters will correspond to 
 a simple dilution process , and the corresponding plot of xm versus x2 will be linear . an example plot ( created using real data ) , which illustrates the method , 
 is shown in figure 1 . due to the linearity the line is then 
 given by13where xm * is the molar property of pure 
 1 ( always water here ) and a is the slope . after taking 
 the partial derivative of eq 13 with respect 
 to n2 , multiplying both sides by n , and taking the limit x2 0 , one finds that14thus , the 
 linear regression of xm provides a and subsequently x2 , the thermodynamically 
 rigorous infinitely dilute 
 partial molar property of the protein plus its counterions . rearrangement 
 of eq 13 , after substituting eq 14 in for a , indicates that this method is 
 equivalent to the apparent molar approach often used experimentally,15although we have seen a 
 few molecular dynamics 
 studies that use this approach to calculate v2 , to our knowledge , it has never been presented generally or applied 
 for the determination of higher derivatives . illustration of the apparent 
 molar approach used to calculate infinitely 
 dilute partial molar properties . a protein is simulated in various 
 box sizes so that a plot of xm vs x2 may be achieved from which x2 may 
 be calculated . ( simulation boxes not to scale . ) in using eq 13 , one has a choice of 
 fitting 
 both xm * and a , or of fixing the y - intercept and only fitting the slope . we have chosen to fix the y - intercept for the following two reasons : ( 1 ) so that we 
 could use the same values of xm * regardless of the protein under 
 study and ( 2 ) because we have performed the pure water simulation 
 for four times longer than the simulations of systems containing proteins , 
 in an effort to increase the precision on the molar properties of 
 water . it should be noted that one could use this method with any 
 number of components as long as the ratio of 1:3: ... :n is held constant and only the ratio of 2:(n the advantages of the above approach include the 
 following:(i)it involves an interpolation of data 
 between pure water and some value of x2 , not an extrapolation.(ii)it is thermodynamically exact.(iii)no subjective definitions are required.(iv)it can be applied to extract any 
 partial molar property for any solute in any solvent system regardless 
 of the number of components.(v)partial molar properties can be assigned 
 to specific conformations . it involves an interpolation of data 
 between pure water and some value of x2 , not an extrapolation . it can be applied to extract any 
 partial molar property for any solute in any solvent system regardless 
 of the number of components . however , regardless 
 of the quality of the approach , in computer 
 simulations , the ultimate arbiter remains the quality of the force 
 field and the need for sufficient sampling of the relevant phase space . specifically , 
 there will be errors in our results due to a classical treatment of 
 those x which have significant quantum mechanical 
 contributions , such as the constant pressure heat capacity . the disadvantage of the approach is that , as alluded to above , 
 multiple simulations are required to obtain a good fit for the calculation 
 of x2. since it is a linear fit , in principle , 
 only two compositions are required , namely , a pure water box and one 
 protein system . however , if any noise is present in the results , it 
 can be better identified when multiple compositions are used . furthermore , 
 systematic deviations from linearity in the plot of xm versus x2 would indicate 
 that a simulation box was too small , such that bulk water was not 
 present . we illustrate the above 
 method using basic pancreatic trypsin inhibitor 
 ( 58 amino acids , pdb i d 5pti ) , abbreviated as bpti , and hen egg white lysozyme 
 ( 129 amino acids , pdb i d 4lzt ) , abbreviated as hew lysozyme . these proteins are 
 relatively small , which allows smaller simulation boxes to be used , 
 and experimental data is available for several of the partial molar 
 properties for comparison with the simulated results . the initial 
 structures for bpti and hew lysozyme were each centered in a series 
 of truncated octahedron boxes where the distances between any two 
 parallel box faces were 8 , 9 , 10 , 11 , and 12 nm . the titratable residues 
 were protonated to correspond to a neutral ph and counterions were 
 added to ensure neutral systems ( 6 cl for bpti 
 systems , 9 cl for hew lysozyme systems ) . a pure 
 water system , consisting of a truncated octahedron box ( with the distance 
 between any two parallel box faces equal to 7 nm ) all 
 simulations were performed using classical md techniques at 
 a temperature of 300 k and pressure of 1 bar unless otherwise noted . 
 in an effort to ensure that the simulations would sample from the npt ensemble and the corresponding system fluctuations would 
 be correct , the nos buff 
 force field ( http://kbff.chem.k-state.edu , kbff ) in explicit solvent with the spc / e water 
 model . following 25 ns of production simulation for each system size , the average c root mean squared deviation ( rmsd ) 
 for the final bpti ( lysozyme ) structures , after a translational and 
 rotational fit to the initial structure , was 0.22 nm ( 0.27 nm ) . these 
 values correspond to the average over the final structure from each 
 of the different system sizes . due to the promising initial results for bpti and hew lysozyme , 
 and due to our overarching questions , we also asked whether the method 
 could be used to distinguish between different conformations of a 
 protein , that is , if it could be used to determine ndx2. 
 to test this , we simulated a native and an extended denatured conformation 
 of trp - cage ( pdb i d 2jof , construct tcb10 ) using the kbff models ( see the supporting information ) . the system setup was the same as for bpti and hew lysozyme , 
 except that the distance between any two parallel box faces was 7 , 
 8 , 9 , and 10 nm for the native conformation , and 8 , 9 , and 10 nm for 
 the denatured conformation . the simulation protocol for trp - cage was 
 similar to that for bpti and hew lysozyme . the specific trp - cage construct 
 used here is neutral , so no counterions were added to the systems . 
 following 25 ns of production simulation , the average c rmsd for the final trp - cage native structures , after a translational 
 and rotational fit to the initial structure , was 0.17 nm . in 
 addition to modeling trp - cage using the kbff in explicit solvent 
 with the spc / e water model , we performed identical trp - cage simulations 
 with the amber99sb force field and the tip3p water model to allow 
 for a force field comparison . a tip3p 
 truncated octahedron box , with the distance between any two parallel 
 box faces equal to 7 nm , was used to calculate the properties of pure 
 water for the amber99sb simulations . following 25 ns of production 
 simulation , the average c rmsd for the final trp - cage 
 native structures , after a translational and rotational fit to the 
 initial structure , was 0.09 nm using amber99sb . the denatured 
 conformation was obtained after clustering a 100 ns , 500 k , nvt simulation using the amber99sb force field . subsequent 
 simulations of this structure using kbff and amber99sb resulted in 
 the denatured conformations shown . the three structures for each force 
 field correspond to the structures from the three compositions ( i.e. , 
 the three system sizes ) studied . the initial and final denatured conformations are shown in 
 figure 2 . since the experimental ndx2 corresponds 
 to a difference between properties of an ensemble of native and denatured 
 states , our goal here was not necessarily to match the experimental 
 value , although certainly that is of long - term interest . in this study , 
 we instead specifically sought to establish how precisely the differences 
 could be calculated , and to see if we achieved the correct sign and 
 order of magnitude for the values using only two conformations . the properties um , km , hm , and vm were extracted from 
 the gromacs energy file since they simply correspond to time averages 
 of system properties . the other properties ( cp , m , t , 
 and p ) were calculated from their 
 bulk system fluctuations definitions in eq 7 . the standard deviation for each property was calculated using block 
 averages where the block size was 5 ns . the system fluctuations 
 displayed a large standard deviation ( see results 
 and discussion ) . in an effort to increase the signal - to - noise 
 ratio , we additionally calculated these properties using the expressions 
 in eq 6 ( where the standard state superscripts 
 should be removed ) . to do this , each system was run at four additional 
 temperatures ( 290 , 295 , 305 , and 310 k ) and four additional pressures 
 ( 250 , 500 , 750 , and 1000 bar ) for 10 ns per new state point . the cp , m , t , and p were then obtained 
 by fitting the simulated molar volumes or enthalpies to a polynomial 
 function and taking the derivatives of the polynomial at the temperature 
 or pressure of interest ( 300 k or 1 bar ) . the standard 
 deviations for each property using this  polynomial fitting 
 method  were calculated using block averages where the block 
 size was again 5 ns . although the production simulation lengths for 
 the system fluctuation approach ( 25 ns ) and polynomial fitting approach 
 ( 10 ns per state point ) are not equal and correspond to 25 ns per 
 property for the system fluctuation approach and 50 ns per property 
 for the polynomial fitting approach , the reported error estimates 
 ( see results and discussion ) indicate that 
 the polynomial fitting method greatly reduced the noise beyond the 
 type of noise reduction one would see by doubling the simulation length 
 alone . figure 3 shows that ensemble average properties 
 were well behaved for both bpti and hew lysozyme . in contrast , there 
 was significant noise in the three properties that are obtained from 
 fluctuating properties . this is clearly illustrated by the scattering 
 of the filled points around the solid lines in figure 4 . it is well - known that the isothermal compressibility 
 of proteins 
 is less than that of pure water . lines 
 with negative slopes in the top panel of figure 4 correspond to protein partial molar isothermal compressibilities 
 that are less than the isothermal compressibility of pure water . if 
 the slope is too negative ( see , e.g. , the system fluctuation result 
 for bpti ) , then the compressibility of the protein will actually become 
 negative ( table 2 ) . system properties used 
 to calculate the infinitely dilute partial 
 molar properties of basic pancreatic trypsin inhibitor ( bpti ) and 
 hen egg white ( hew ) lysozyme in pure water at 300 k and 1 bar according 
 to eq 14 . circles : 
 the property at each composition . lines : linear fits ( fixed y - intercept ) through the set of points . the properties shown 
 correspond to time averages from the simulations : the molar internal 
 energy , um ( kjmol ) , 
 kinetic energy , km ( kjmol ) , volume , vm ( nm  
 10 ) , and enthalpy , hm ( kjmol ) . system properties used 
 to calculate the infinitely dilute partial 
 molar properties of basic pancreatic trypsin inhibitor ( bpti ) and 
 hen egg white ( hew ) lysozyme in pure water at 300 k and 1 bar according 
 to eq 14 . results 
 from the bulk system fluctuations are shown as solid lines with filled 
 circles , while results from the polynomial fitting method are shown 
 as dotted lines with open circles . the property at each composition 
 is shown as a circle , and the linear fit through the set of points , 
 in which the y - intercept was fixed , is shown as a 
 line . the top panel corresponds to vmt ( 10 nm bar ) , the middle panel to cp , m ( j mol k ) , and the bottom panel to vmp ( 10 nm k ) . error estimates are reported in 
 parentheses and correspond to the standard deviation from 5 ns block 
 averages . the experimental hew lysozyme volume was measured 
 using a solution of isoionic lysozyme in distilled and deionized water . the volume of a single spc / e water 
 molecule is 0.03 nm , and 1 nm molecule = 602.2 cm mol . error estimates are reported in parentheses and correspond 
 to the standard deviation from 5 ns block averages . experimental bpti 
 heat capacity value from makhatadze et al . , and thermal expansion value from lin et al . experimental lysozyme heat capacity values from gekko and noguchi and yang and rupley , isothermal compressibility value from gekko and hasegawa , and thermal expansion value from gekko and noguchi . tables 1 and 2 show the values and standard deviations of the x2 obtained 
 for bpti and hew lysozyme . while the internal 
 energy , kinetic energy , and enthalpy of the proteins are not comparable 
 to experiments , it may be of interest to study trends in these properties 
 with , for example , protein mass , for a given force field . one application 
 of such a study might be that these trends could provide an additional 
 metric for the characterization of protein force fields . it is possible 
 that such a characterization , for a given force field , could be used 
 to distinguish between folded and unfolded conformations using thermodynamic , 
 rather than structural , criteria . the simulated protein kinetic energy 
 agreed with the predicted equipartition result in every case , calculated 
 as k2 = 1/2kbt(3natoms  nbonds ) , where natoms is the number of protein + ion atoms and nbonds is the number of constrained bonds in the protein . the simulated bpti volume , 8.164(8 ) nm , was larger 
 than the approximate value calculated based upon the amino acid composition , 
 7.8 nm . this difference of 
 0.364 nm corresponds to approximately 12 waters ( the volume 
 of an spc / e water molecule is 0.03 nm ) . the simulated 
 hew lysozyme volume , 18.00(3 ) nm , was larger than the 
 experimental volume , 16.92(2 ) nm ( experimental measurement 
 was of isoionic lysozyme in completely deionized and distilled water 
 at 25 c ) . this difference of 1.08 
 nm for hew lysozyme corresponds to approximately 36 waters . 
 for both bpti and hew lysozyme , the disagreement with the literature 
 results is partly due to the inclusion of counterion volumes in our 
 simulated values ( 6 for bpti , 9 for hew lysozyme ) . for hew lysozyme , 
 this is clearly not the only issue . using the system fluctuation 
 method , the fluctuating properties 
 for both proteins displayed too much noise to determine if there was 
 good agreement with experiment or not ( table 2 ) . to increase the signal - to - noise ratio , we used a second approach , 
 in which we looked at the system volume and enthalpy as a function 
 of pressure and or temperature and calculated the t , cp , m , and p using their partial derivative 
 definitions instead of their fluctuation definitions ( see theory ) . this greatly reduced the error ( shown visually 
 in figure 4 and tabulated in table 2 ) . since our simulations 
 were classical and there are significant quantum corrections to the 
 heat capacity , perfect agreement with experiment , even if achieved , 
 may have been for the wrong reasons . specifically , real molecules 
 store energy in low frequency bonds , whereas all bonds were constrained 
 in our simulations . second , in real molecules , only the low frequency 
 bond angles would store energy , whereas all angles ( except those in 
 water ) are flexible and store energy in our simulations . quantum corrections 
 could be made , but they are difficult 
 for proteins and have not been made here . the thermal expansion 
 and isothermal compressibility were in reasonable 
 agreement with experiment , but displayed large error bars . the magnitude 
 of the errors on the fluctuating properties is a downside of our approach 
 that has not been observed using more subjective methods . literature values are available for all of the pure water 
 ( spc / e ) 
 properties . however , since it is critical that all of the simulation 
 conditions are consistent , we simulated a pure spc / e system instead 
 of simply quoting previously reported spc / e properties . the reason 
 for the difference between the pure spc / e values obtained using the 
 system fluctuation method versus the polynomial fitting method for 
 the isothermal compressibility is currently unknown ( figure 4 ) . however , the value obtained using the polynomial 
 fitting method is in better agreement with experiment [ polynomial 
 fitting , 4.48(2 )  10 bar ; system fluctuation , 4.75(3 )  10 bar ; experiment , 4.525  
 10 bar ] . the results 
 for the trp - cage simulations with both force fields 
 and with both methods ( the system fluctuation and the polynomial fitting 
 method ) are summarized in tables 3 and 4 . as shown in table 3 , the 
 kbff results suggest that the native state has a larger volume than 
 the denatured state by a value of approximately one water molecule 
 or 1% of the protein volume , regardless of whether or not 
 the denatured conformation was position restrained . this is the typical 
 sign observed for the volume change of a protein upon denaturation . error estimates are reported in 
 parentheses and correspond to the standard deviation from 5 ns block 
 averages . the kinetic energy contributions agreed with the equipartition 
 results and were the same , within statistical uncertainty , for the 
 native and denatured confirmations . error estimates 
 are reported in 
 parentheses and correspond to the standard deviation from 5 ns block 
 averages , excluding the 1 s simulation , which corresponds to 
 the standard deviation on 333 ns block averages . for the 1 
 s simulations , the error on the y - intercept 
 ( pure water value ) was taken as zero . sf : system fluctuations , pf : 
 polynomial fitting . the 
 kbff simulation produced the correct sign for the change in 
 enthalpy upon denaturation . the enthalpy of the denatured state was 
 more positive , or unfavorable , than the enthalpy of the native state , 
 as expected for heat denaturation . the average 
 kbff result was in better agreement with experiment when the position 
 restraints were removed , although the difference was not statistically 
 significant and the native state was slightly overstabilized in both 
 cases . since the starting conformation for both simulations was taken 
 from a 500 k simulation using the amber99sb force field ( see the supporting information ) , better agreement with 
 experiment may be achieved if the kbff simulations were extended to 
 allow for more conformational rearrangements . the amber99sb result 
 was within the error of the experimental value with or without position 
 restraints on the denatured conformation . globular protein heat 
 capacity changes upon denaturation are typically 
 positive ; a positive ( negative ) change is considered to be due to 
 apolar ( polar ) solvation upon denaturation . since trp - cage is considered globular and its name stems from its 
 buried tryptophan side chain , one might predict that the sign of ndcp,2 would be positive . indeed , the experimental ndcp,2 for the original trp - cage construct ( tc5b ) was reported 
 to be small and positive ( 0.3  0.1 kj mol k ) . were small and negative ( ndcp,2 = 0.2  0.1 kj mol k ) for the trp - cage construct simulated here 
 ( tc10b ) , indicating that considering 
 tryptophan exposure is not enough to predict the correct sign . while the statistical uncertainty in our results was too great 
 to determine whether or not the correct sign was obtained , our results 
 may suggest that the kbff and amber99sb simulations provide opposite 
 signs ( negative and positive , respectively ) . it is encouraging that 
 the kbff simulations may be able to detect this atypical ndcp,2 sign . the positive sign for the amber results reported 
 here agrees with the recent 1 s per replica remd simulations 
 of the tc10b construct by english and garcia , who also used the amber99sb / tip3p 
 force fields and reported ndcp,2 = 0.35 
  0.51 kj mol k. the results for the isothermal compressibility 
 and thermal expansion 
 changes upon denaturation are too noisy to determine the simulated 
 sign ; however , the error was greatly reduced using the polynomial 
 fitting method . visualization of the denatured trp - cage trajectories 
 showed that 
 the kbff and amber99sb denatured simulations without position restraints 
 behaved differently from each other . the kbff denatured trp - cage conformation 
 remained more fully extended throughout the 25 ns , whereas the amber99sb 
 denatured trp - cage collapsed quickly and remained collapsed for most 
 of the simulation ( although it became more extended at the very end 
 of the simulation ) . the final snapshot of the denatured simulations 
 is shown in figure 2 . despite this , the two 
 force fields gave similar enthalpy differences [ experiment , ndh2 = 
 65(2 ) kj mol ; kbff , ndh2 = 90(20 ) 
 kj mol ; amber99sb , ndh2 = 50(20 ) 
 kj mol ] . certainly , the simulations did not access 
 all the denatured conformations with only 25 ns of production , so 
 it was remarkable that the ndh2 values were so close to 
 the experimental value . to test whether the different denatured 
 conformations that were 
 sampled in the simulations of trp - cage with these two different force 
 fields accounted for the different enthalpy differences , we ran a 
 new set of denatured simulations for both ffs using soft harmonic 
 position restraints ( 100 kj mol nm ) on the c atoms of trp - cage . the average position 
 restraint energy was small , approximately 2 j mol , so we did not repeat the simulations of the native trp - cage with 
 position restraints . thus , we assumed that any artifacts from the 
 addition of the position restraints were negligible when calculating 
 the ndx2 using simulations of the native structure 
 without position restraints . the ndx2 results 
 were the same with and without position restraints to within the statistical 
 error . this is consistent with a two - state folding pathway in which 
 the energy landscape is relatively flat for the ( high free energy ) 
 unfolded conformations and drops down into one minimum for the folded 
 conformation . remarkably , the enthalpy differences remained essentially 
 correct when calculated from the difference between the native conformation 
 and just one restrained denatured conformation . references therein have 
 reported that typical magnitudes for the volume change associated 
 with protein unfolding are < 0.5% of the protein volume . using pressure perturbation calorimetry , lin 
 determined the ndv2 value was 0.06% of the volume 
 of native bpti in water at ph 4 and 0.08% of the volume of 
 native hew lysozyme in water at ph 2.5 . the errors on our ndv2 value for bpti and hew 
 lysozyme were 0.1% and 0.2% ; thus , simulation times longer than the 
 currently used times would be necessary to calculate ndv2 with 
 precision for these proteins . the corresponding ndp,2 values reported by lin were 0.5  10 k and 1.4  10 k for bpti and hew lysozyme , respectively . to calculate ndp,2 with sufficient precision using the polynomial fitting method , we 
 would need to reduce the error on p,2 by 
 a factor of 8 for bpti , which would require simulations approximately 
 64 times as long , that is , 640 ns at each temperature , assuming that 
 the error  tsim . with the current simulation time , the p,2 value for hew lysozyme was at the minimal precision necessary for 
 detection of ndp,2 , if simulations of denatured 
 conformations were performed . in a compendium of protein thermodynamic 
 properties makhatadze has reported that bpti ndh2 = 130 kj 
 mol and ndcp,2 = 3.0 
 kj k mol at 298 k. these numbers are larger than the uncertainty 
 on the native enthalpy and native heat capacity calculated here , 10 
 kj mol and 0.4 kj mol k , respectively . the same source reports that hew lysozyme 
 ndh2 = 242 kj mol and ndcp,2 = 9.1 kj mol k at 298 k. again , these numbers are 
 larger than the uncertainty on the native hew lysozyme enthalpy and 
 heat capacity determined here , 50 kj mol and 0.09 
 kj mol k , respectively . thus , 
 we predict that similar length simulations of either of these proteins 
 in their denatured conformations would allow for ndh2 and 
 ndcp,2 to be calculated with statistical precision . the challenge in these situations , heightened as the number of residues 
 in the protein of interest increases , would be in choosing protein 
 conformations that were  representative  of the denatured 
 state . although the trp - cage results from these simulations surprisingly 
 showed that the values were not sensitive to the protein conformation , 
 much more work on different proteins would be necessary before confirming 
 that behavior as a general conclusion . . we could have simulated all of the systems for even longer 
 periods 
 of time . however , even a multiplication of the production simulation 
 length by four would only result in a noise reduction by a factor 
 of 2 ( assuming that the error  tsim ) , which would not be sufficient in most cases 
 to determine these properties with sufficient precision . so , in an 
 effort to estimate how long one would need to simulate trp - cage , in 
 order to obtain the heat capacity and thermal expansion with sufficient 
 precisely , we extended the amber99sb native trp - cage simulations to 
 1 s at temperatures of 290 , 295 , 300 , 305 , and 310 k. since 
 we had already established the linear behavior of xm versus x2 for the 7 , 8 , 
 9 , and 10 nm box sizes , we only extended the simulations for the 8 
 nm box . extension of a simulation from 10 ns to 1 s is an increase 
 of tsim by a factor of 100 ; thus , we would 
 have predicted a reduction in the noise by a factor of 10 . as shown 
 in table 4 , the errors were actually reduced 
 by a factor of 5 for cp,2 and 
 only a factor of 2.5 for p,2. we have investigated the ability of a simple method to extract 
 infinitely dilute partial molar properties of biomolecules from molecular 
 dynamics simulations using basic pancreatic trypsin inhibitor and 
 hen egg white lysozyme as our test case proteins . we then used the 
 method to distinguish between the protein volume and enthalpy of a 
 native and a denatured trp - cage conformation with sufficient precision . the method is most similar to that used recently by smit to determine 
 the enthalpy of a lipid aggregate , and to shing and chung s 
 work in the late 1980s for the infinitely dilute partial molar volume 
 and internal energy of a lennard  jones solute , but to our knowledge it has never been used or presented for any 
 general solute and for any general partial molar property . the strengths 
 of the approach are that it is thermodynamically exact , no ambiguous 
 parameters are introduced , and knowledge of the specific conformations 
 that contributed to the output values is retained . the method is general , 
 and this approach may be used to calculate any partial molar property 
 of any solute in any solution . it could be employed in the future 
 as a test of the quality of a ff . the current limitation of 
 the approach is the high level of noise 
 in the isothermal compressibility and thermal expansion coefficient 
 results , and the moderate level of noise in the isobaric heat capacity 
 results . to increase the signal - to - noise ratio , a polynomial fitting 
 method was used . clearly , the system fluctuation approach was not 
 the method of choice for the calculation of the fluctuating properties , 
 and the polynomial fitting method provided higher precision results . however , we were still unable to precisely distinguish between native 
 and denatured trp - cage conformations heat capacities , thermal 
 expansions , and isothermal compressibilities . based upon our 
 results , we predict that for a protein the size 
 of hew lysozyme ( 129 residues ) , the following amounts of time are 
 needed to calculate specific thermodynamic properties with statistical 
 precision : ndv2 , 200 ns ; ndh2 , 
 achieved necessary precision using 25 ns simulation ; ndcp,2 , achieved necessary precision using 10 ns simulations 
 at each of 5 temperatures ; ndt,2 , difficult 
 to predict due to the large noise ; ndp,2 , achieved the 
 necessary precision using 10 ns simulations at each of 5 temperatures . 
 currently , the methods that rely upon subjective definitions , if they 
 can be trusted , are more competitive than our approach for the isothermal 
 compressibility due to the high amounts of noise in the current method . 
 continued increases in computational power should make the method 
 presented here more usable , with ease , over time . lastly , the 
 method presented here does not provide a means of directly 
 decomposing the thermodynamic properties into contributions , from 
 specific surface groups for example . the solute  solute and 
 solute  solvent energy terms decay as a function of distance 
 away from the solute , which makes them directly calculable . since 
 the solvent  solvent contribution does not decay , it must be 
 solved for indirectly . this is a commonly encountered issue ; see , 
 for example , mccammon and co - workers . we are currently working to develop and test a method to calculate 
 all of the above properties that uses expressions from fluctuation 
 solution theory and is decomposable , analogous to the group contribution 
 decomposition for partial molar volumes already achieved using traditional 
 kirkwood  buff theory .
a detailed 
 understanding of temperature and pressure effects on 
 an infinitely dilute protein s conformational equilibrium requires 
 knowledge of the corresponding infinitely dilute partial molar properties . 
 established molecular dynamics methodologies generally have not provided 
 a way to calculate these properties without either a loss of thermodynamic 
 rigor , the introduction of nonunique parameters , or a loss of information 
 about which solute conformations specifically contributed to the output 
 values . here 
 we implement a simple method that is thermodynamically 
 rigorous and possesses none of the above disadvantages , and we report 
 on the method s feasibility and computational demands . 
 we calculate 
 infinitely dilute partial molar properties for two proteins and attempt 
 to distinguish the thermodynamic differences between a native and 
 a denatured conformation of a designed miniprotein . 
 we conclude that 
 simple ensemble average properties can be calculated with very reasonable 
 amounts of computational power . 
 in contrast , properties corresponding 
 to fluctuating quantities are computationally demanding to calculate 
 precisely , although they can be obtained more easily by following 
 the temperature and/or pressure dependence of the corresponding ensemble 
 averages .