Computational investigation of dynamical transitions in Trp-cage miniprotein powders
Computational investigation of dynamical transitions in Trp-cage miniprotein powders
Sang Beom Kim
Devansh R. Gupta
Pablo G. Debenedetti
OPEN We investigate computationally the dynamical transitions in Trp-cage miniprotein powders, at three levels of hydration: 0.04, 0.26 and 0.4 g water/g protein. We identify two distinct temperatures where transitions in protein dynamics occur. Thermal motions are harmonic and independent of hydration level below Tlow ? 160 K, above which all powders exhibit harmonic behavior but with a different and enhanced temperature dependence. The second onset, which is often referred to as the protein dynamical transition, occurs at a higher temperature TD that decreases as the hydration level increases, and at the lowest hydration level investigated here (0.04 g/g) is absent in the temperature range we studied in this work (T ? 300 K). Protein motions become anharmonic at TD, and their amplitude increases with hydration level. Upon heating above TD, hydrophilic residues experience a pronounced enhancement in the amplitude of their characteristic motions in hydrated powders, whereas it is the hydrophobic residues that experience the more pronounced enhancement in the least hydrated system. The dynamical transition in Trp-cage is a collective phenomenon, with every residue experiencing a transition to anharmonic behavior at the same temperature.
Although proteins are often described with static representations obtained from X-ray crystallographic data, the
dynamics of proteins are essential to their functionality, as is clearly illustrated, for example, in the case of ligand
binding and ion channel regulation1,2. The characteristic amplitude of protein motions, often measured by
atomistic mean-square fluctuations (MSF) of the protein atoms, increases linearly with temperature up to approximately
180?240 K3,4. At this temperature the amplitude exhibits a sharp transition to a non-linear temperature
dependence, and this onset of anharmonicity is referred to as the protein dynamical transition (PDT)3. Experimental
studies report impaired protein functionality at temperatures below the dynamical transition temperature (TD)5?7.
For a given protein, the PDT shifts to higher temperatures with increasing solvent viscosity, which is
sensitive both to hydration level and solvent composition8?12. Despite many experimental and computational
studies4,8?14,14?20, the underlying physical basis of the PDT is still under debate. Some consider the main
contribution to this transition to come from the activation of protein side-chain dynamics due to the increased
translational and rotational dynamics of water14,17,18. It has also been suggested that the PDT is connected to the
fragile-to-strong transition of the hydration water and the hypothesized liquid-liquid transition, and corresponds
to the crossing of the Widom line19,20.
It has been found that there exists another transition in protein dynamics at a temperature lower than TD4,13?16.
Protein dynamics exhibit an onset of enhanced motion at this transition temperature, which we refer to as Tlow in
this paper. The Tlow is known to be dependent on the type of protein, typically ranging from 100 K to 180 K21?23.
The activation of methyl group motions21,22 and proline puckering transitions23 have been suggested as the
underlying causes for this transition, but as with the PDT, the physical origin of this transition remains unclear.
In order to obtain atomic-level physical insight into these low-temperature transitions, we present a
simulation study of protein powder systems with varying degrees of hydration level. Molecular dynamics (MD)
simulations provide the appropriate spatial and temporal resolution to probe protein dynamics on a microscopic scale,
complementing experimental studies. As a model protein, we used a 20-residue miniprotein Trp-cage, which is
one of the smallest synthetic peptides that show protein-like secondary and tertiary structures (PDB ID: 1L2Y)24.
At ambient conditions, Trp-cage has a well-defined hydrophobic core with both ?-helix and 310-helix
structures24. The small size of Trp-cage makes it an ideal candidate for our study, where the simulations of multiple
Trp-cages in a unit cell are needed to model powder-like environments. We chose to perform simulations of such
an environment, rather than of a system composed of a single protein unit in solution because the powder system
provides more information on protein dynamics directly relevant to practical applications, such as solid-state
The powder systems, each comprised of 16 Trp-cages, were prepared at three different hydration levels. We
refer to these systems as P?0.40, P?0.26, and P?0.04 in this paper, where the number denotes the hydration level
in g water/g protein (g/g). The three hydration levels were chosen to represent fully hydrated, partially hydrated,
and dehydrated powders. Trp-cage is fully hydrated at a hydration level of 0.40 g/g, according to its water sorption
isotherm26, and 0.04 g/g represents the approximate amount of residual, strongly bound water present in typical
freeze-dried protein powders27. We systematically identify the two transitions and their dependence on the
hydration level. We also show the effect of local hydrophobicity on residue-level protein dynamics.
Figure?1 shows the average mean-square fluctuation (MSF) of protein heavy (non-hydrogen) atoms as a function
of temperature, for each hydration level. Below ?150 K, the protein dynamics are essentially harmonic,
exhibiting a linear increase in the average MSF with temperature. In this range of temperatures, the magnitudes of the
average MSFs are essentially identical among all hydration levels considered here. Thus, the average MSFs of the
protein heavy atoms are not affected by the level of hydration at temperatures below Tlow.
Figure?2a shows the average MSFs in the temperature range T = 100?200 K, with linear fits to the 100?150 K
regions. At T ? 160 K, the average MSF starts deviating from the low-temperature linear increase with
temperature. This temperature corresponds to the first transition temperature (Tlow), and the powders with varying
hydration levels have a common value of Tlow, within the numerical accuracy of our calculations. The independence of
Tlow on hydration level is consistent with previous studies21,28.
Figure?2b shows the average MSF in semi-logarithmic scale as a function of temperature. The two solid lines
represent linear fits to the data below Tlow (?160 K) and between Tlow and 210 K. This indicates that the Trp-cage
powders exhibit a linear T-dependence of the average MSF above Tlow as well, but with an increased slope. With
the exception of P?0.04 system that shows no PDT up to 300 K, the average MSF then starts to increase
nonlinearly at T ? 220?240 K. We refer to this temperature as the TD. Above TD, the amplitude of protein motions is
significantly enhanced with increasing hydration level.
It can be nontrivial to determine TD by locating the temperature at which the MSF starts increasing
nonlinearly, especially for less-hydrated systems. This is because the onset of the anharmonic increase in dynamics
becomes less pronounced as the hydration level decreases. It has been suggested that the PDT is related to the
glass transition of the partially hydrated protein system29. Accordingly, we investigate the temperature
dependence of the enthalpy, which should, according to this hypothesis, exhibit a discontinuous change in slope at TD,
corresponding to a jump in the heat capacity30. Figure?3 shows the enthalpy of each powder system as a function
of temperature. With the exception of the dehydrated system (P?0.04), two linear regimes with different slopes
are indeed present. The question of whether this change of slope corresponds to a true glass transition is one that
we do not address here; we simply note that this particular aspect of glassy behavior30 is present and, remarkably,
we find these effective glass transition temperatures to be in very good agreement with the TD found in Fig.?2b
(dashed lines). TD shifts to higher temperatures as the protein becomes dehydrated, confirming for the Trp-cage
powders investigated here the experimentally observed dependence of the PDT upon the hydration level8?12. The
marginally hydrated, P?0.04 system, exhibits a single linear regime for the enthalpy as a function of T, due to the
lack of the PDT up to the highest temperature investigated here, T = 300 K. Finally, we note that the heat capacity
jump, ?Cp, increases with the hydration level.
In order to understand how each residue contributes to the overall dynamics of the Trp-cage powders, we
show the MSF of each residue in Fig.?4. It can be seen that certain residues exhibit especially high MSF values
compared to the rest. The termini residues (N1 and S20) have high MSF values because they are only bonded to
one other residue and thus experience comparatively less restraints. In the hydrated powders (i.e. P?0.40 and
P?0.26), the non-terminus residues with high MSF values (Q5, K8, D9, S13 and R16) have the common trait of
being hydrophilic, due to either nonzero charges (K8, D9, and R16) or polar functional groups that can form
hydrogen bonds with water (Q5 and S13).
The dehydrated system (i.e. P?0.04), however, exhibits opposite behavior. Figure?5a shows the Trp-cage unit,
with each residue colored according to the hydrophobicity scale of Eisenberg et al.31. Figure?5b compares the
powders with hydration levels 0.4 and 0.04 g/g by coloring each residue according to the magnitude of its MSF.
As the temperature is increased, the MSF of the hydrophilic residues increases relative to that of the hydrophobic
residues in the hydrated powder. In contrast, it is the hydrophobic residues (e.g. L2, W6 and L7) of the dehydrated
powder that experience enhanced MSF upon heating. This enhancement of fluctuations in the hydrophobic
We have investigated the dynamical transitions of Trp-cage powders with varying hydration levels (h = 0.40, 0.26,
and 0.04 g/g). We identified two temperatures where transitions in the temperature-dependent protein dynamics
occur. The average MSF of the protein heavy atoms increases linearly with temperature at low temperatures (100?
150 K). The first transition at Tlow (?160 K) is characterized by a sudden change in the slope of this harmonic
(linear) behavior. The second transition, PDT, occurs at a higher temperature (TD) where the MSF starts to depend
nonlinearly upon temperature. We showed that the TD can be precisely determined by locating the ?calorimetric
glass transition temperature? of the protein/water system, where a sudden jump in the heat capacity occurs. We
find this heat-capacity jump to be more pronounced the larger the hydration level. The Tlow and the MSF below
Tlow are identical for powder systems of different hydration levels. In contrast, dehydration shifts the TD to a
higher temperature, while greatly suppressing the protein dynamics at T > TD.
As the powder systems are heated, the dynamics of methyl-group hydrogens, as measured by the MSF,
increases exponentially with T up to Tlow, and linearly for T > Tlow, without any signs of a transition at TD. In
contrast, the MSF of non-methyl hydrogens increases linearly up to TD and nonlinearly for T > TD. This linear
increase in the MSF of non-methyl hydrogens below TD occurs without any slope change at Tlow. We have
therefore identified distinct contributions from the methyl and non-methyl groups to the transitions at Tlow and TD.
Our findings on the temperature dependence of the methyl-group dynamics of the Trp-cage, however, are in
contrast to previous studies21,22 where the activation of the methyl groups was suggested as the cause of the transition
at Tlow. This suggests that important aspects of low-temperature protein dynamics may resist generalization across
different individual proteins.
All 20 residues of the Trp-cage display the two transitions in their dynamics at the same Tlow and TD as the
whole protein. The magnitude of their MSFs, however, depends largely on the degree of
hydrophobicity/hydrophilicity of each residue. As the hydrated system (P?0.40) is heated, the characteristic motions of the hydrophilic
residues are enhanced relative to those of the hydrophobic residues. The opposite is true for the dehydrated
system (P?0.04) where the hydrophobic residues exhibit more pronounced increases in their MSFs. This points to
hydration as an important process variable in modulating the temperature stability of solid-state pharmaceutical
We have presented a computational study of dynamical transitions in Trp-cage powders with varying
hydration levels. It will be interesting to perform similar simulation studies with larger proteins or other biomolecules,
such as RNAs and DNAs, in order to assess the generality of our findings. Another potentially fruitful avenue of
inquiry will be to study a more complex protein matrix that consists not only of protein and water but includes
also cosolutes, such as carbohydrates, that are commonly present in pharmaceutical formulations. By
investigating the effects of sugar molecules on the residue-level dynamics of the simulated proteins, such studies would
complement previous investigations on the shift of the dynamical transition temperatures due to the presence of
System Preparation. The powder system at a hydration level of 0.40g/g was prepared by following the
procedures described in detail in ref. 26. It contains 16 Trp-cages that are randomly translated and rotated, 771 water
molecules, and 16 chloride ions. The net positive (+1e) charge of each Trp-cage at neutral pH was balanced by
the negatively charged chloride ion (?1e). Powders at lower hydration levels were prepared by dehydrating this
powder, through cycles comprising the removal of one water molecule and the relaxation of the resulting system
through 200 ps of NPT MD simulation at 300 K and 1 bar26. The water molecule to be removed was randomly
chosen and accepted/rejected by standard Metropolis criteria based on the associated Boltzmann factor (e???U),
where ? is the inverse of the product of Boltzmann?s constant and temperature, and ?U is the change in
configurational energy that would result from removing the water molecule in question. In this way, water molecules that
are more strongly bound to protein units are less likely to be removed, preventing a potentially large perturbation
to the protein structure.
Molecular Dynamics Simulation. The GROMACS33?36 package was used for all MD simulations. The
leap-frog algorithm was used to integrate the equations of motion, with a time step of 1 fs. Temperature and
pressure were controlled using Nos?-Hoover thermostat37,38 with a 0.1 ps time constant and the Parrinello-Rahman
barostat39,40 with a 1 ps time constant, respectively. Anisotropic pressure coupling ensured that fluctuations in the
dimensions of the orthorhombic simulation box were independent of each other. Periodic boundary conditions
were applied in all three dimensions. We truncated the short-range interactions at 1 nm and applied the
standard long-range dispersion corrections for the energy and pressure41. The reciprocal part of the Ewald sum for
long-range electrostatics was calculated using the smooth particle mesh Ewald method42. The linear constraint
solver algorithm (LINCS)43,44 and SETTLE45 were used to constrain all bonds in the protein and water molecules,
respectively. Proteins and water molecules were modeled using Amber ff03w46,47 and TIP4P/200548 force fields,
The systems were equilibrated at 300K and 1 bar initially, through 5 ns of NVT MD, followed by 5 ns of NPT
MD. Each prepared system was replicated and equilibrated to low temperatures (100?300 K) by performing an
NPT MD simulation during which the temperature was decreased linearly from 300 K to the temperature of
interest at a rate of 4 K/ns. Variations of this cooling rate (1, 4, and 10 K/ns) were tested for one of the powders
(h = 0.40 g/g), and identical protein dynamics were observed. We then performed 200 ns of NPT MD and used
the last 150 ns of the trajectory for analysis. We performed the block-averaging analysis49,50 by dividing each
trajectory into 5 blocks, in order to estimate the standard error for each observable we analyzed.
Mean-Square Fluctuation. In order to compute the mean-square fluctuation (MSF), the structure of each
protein (using protein heavy atoms) at all time steps was first aligned to the structure of the corresponding protein
at time = 0. The MSF of each protein heavy atom, i, was then calculated by computing the variance in its atomic
positions from the average position:
1 N F 2
? ???xi(t j) ? xi ??? ,
N F j=1
where NF is the number of frames in a trajectory, xi is the position of an atom i, and ?? is the ensemble average.
The atomic MSF values were then mass-averaged to compute the MSF for each Trp-cage, and the MSF of the
individual Trp-cages were then averaged to yield the average MSF at each temperature.
P.G.D. gratefully acknowledges financial support from the National Science Foundation (Grant No.
CBET1263565). The computations were performed at the Terascale Infrastructure for Groundbreaking Research in
Engineering and Science (TIGRESS), at Princeton University. This work also used the Extreme Science and
Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation (Grant No.
S.K. and P.G.D. designed the research. S.K. performed all the simulations. D.R.G. contributed to some of the early
exploratory simulations. S.K. and P.G.D. analyzed and discussed the results. S.K. and P.G.D. wrote and edited the
Supplementary information accompanies this paper at http://www.nature.com/srep
Competing financial interests: The authors declare no competing financial interests.
How to cite this article: Kim, S. B. et al. Computational investigation of dynamical transitions in Trp-cage
miniprotein powders. Sci. Rep. 6, 25612; doi: 10.1038/srep25612 (2016).
This work is licensed under a Creative Commons Attribution 4.0 International License. The images
or other third party material in this article are included in the article?s Creative Commons license,
unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license,
users will need to obtain permission from the license holder to reproduce the material. To view a copy of this
license, visit http://creativecommons.org/licenses/by/4.0/
1. Teilum , K. , Olsen , J. G. & Kragelund , B. B. Protein stability, flexibility and function . Biochim. Biophys. Acta, Proteins Proteomics 1814 , 969 - 976 ( 2011 ).
2. Allen , T. W. , Andersen , O. & Roux , B. On the importance of atomic fluctuations, protein flexibility, and solvent in ion permeation . J. Gen. Physiol . 124 , 679 - 690 ( 2004 ).
3. Doster , W. , Cusack , S. & Petry , W. Dynamical transition of myoglobin revealed by inelastic neutron scattering . Nature 337 , 754 - 756 ( 1989 ).
4. Zaccai , G. How soft is a protein? a protein dynamics force constant measured by neutron scattering . Science 288 , 1604 - 1607 ( 2000 ).
5. Parak , F. et al. Evidence for a correlation between the photoinduced electron transfer and dynamic properties of the chromatophore membranes from rhodospirillum rubrum . FEBS Lett . 117 , 368 - 372 ( 1980 ).
6. Ferrand , M. , Dianoux , A. J. , Petry , W. & Zaccai , G. Thermal motions and function of bacteriorhodopsin in purple membranes: Effects of temperature and hydration studied by neutron scattering . Proc. Natl. Acad. Sci. USA 90 , 9668 - 9672 ( 1993 ).
7. Rasmussen , B. F. , Stock , A. M. , Ringe , D. & Petsko , G. A. Crystalline ribonuclease a loses function below the dynamical transition at 220 K . Nature 357 , 423 - 424 ( 1992 ).
8. Cordone , L. , Ferrand , M. , Vitrano , E. & Zaccai , G. Harmonic behavior of trehalose-coated carbon-monoxy-myoglobin at high temperature . Biophys. J 76 , 1043 - 1047 ( 1999 ).
9. Fitter , J. The temperature dependence of internal molecular motions in hydrated and dry ?-amylase: The role of hydration water in the dynamical transition of proteins . Biophys. J 76 , 1034 - 1042 ( 1999 ).
10. Tsai , A. M. , Neumann , D. A. & Bell , L. N. Molecular dynamics of solid-state lysozyme as affected by glycerol and water: A neutron scattering study . Biophys. J 79 , 2728 - 2732 ( 2000 ).
11. Paciaroni , A. , Cinelli , S. & Onori , G. Effect of the environment on the protein dynamical transition: A neutron scattering study . Biophys. J 83 , 1157 - 1164 ( 2002 ).
12. Zhou , P. & Labuza , T. P. Effect of water content on glass transition and protein aggregation of whey protein powders during shortterm storage . Food Biophys . 2 , 108 - 116 ( 2007 ).
13. Roh , J. H. et al. Onsets of anharmonicity in protein dynamics . Phys. Rev. Lett . 95 , 038101 ( 2005 ).
14. Zanotti , J.-M. , Gibrat , G. & Bellissent-Funel , M.-C. Hydration water rotational motion as a source of configurational entropy driving protein dynamics . crossovers at 150 and 220 K. Phys. Chem. Chem. Phys . 10 , 4865 - 4870 ( 2008 ).
15. Mazza , M. G. et al. More than one dynamic crossover in protein hydration water . Proc. Natl. Acad. Sci. USA 108 , 19873 - 19878 ( 2011 ).
16. Vural , D. , Hong , L. , Smith , J. C. & Glyde , H. R. Long-time mean-square displacements in proteins . Phys. Rev. E 88 , 052706 ( 2013 ).
17. Tournier , A. L. , Xu , J. & Smith , J. C. Translational hydration water dynamics drives the protein glass transition . Biophys. J. 85 , 1871 - 1875 ( 2003 ).
18. Wood , K. et al. Coincidence of dynamical transitions in a soluble protein and its hydration water: direct measurements by neutron scattering and MD simulations . J. Am. Chem. Soc . 130 , 4586 - 4587 ( 2008 ).
19. Chen , S.-H. et al. Observation of fragile-to-strong dynamic crossover in protein hydration water . Proc. Natl. Acad. Sci. USA 103 , 9012 - 9016 ( 2006 ).
20. Kumar , P. et al. Glass transition in biomolecules and the liquid-liquid critical point of water . Phys. Rev. Lett . 97 , 177802 ( 2006 ).
21. Roh , J. H. et al. Influence of hydration on the dynamics of lysozyme . Biophys. J 91 , 2573 - 2588 ( 2006 ).
22. Wood , K. et al. The low-temperature inflection observed in neutron scattering measurements of proteins is due to methyl rotation: Direct evidence using isotope labeling and molecular dynamics simulations . J. Am. Chem. Soc . 132 , 4990 - 4991 ( 2010 ).
23. Miao , Y. et al. Temperature-dependent dynamical transitions of different classes of amino acid residue in a globular protein . J. Am. Chem. Soc . 134 , 19576 - 19579 ( 2012 ).
24. Neidigh , J. W. , Fesinmeyer , R. M. & Andersen , N. H. Designing a 20-residue protein . Nat. Struct. Biol . 9 , 425 - 430 ( 2002 ).
25. Tarek , M. & Tobias , D. J. The dynamics of protein hydration water: A quantitative comparison of molecular dynamics simulations and neutron-scattering experiments . Biophys. J 79 , 3244 - 3257 ( 2000 ).
26. Kim , S. , Palmer , J. C. & Debenedetti , P. G. A computational study of the effect of matrix structural order on water sorption by trpcage miniproteins . J. Phys. Chem. B 119 , 1847 - 1856 ( 2015 ).
27. Wang , W. Lyophilization and development of solid protein pharmaceuticals . Int. J. Pharm . 203 , 1 - 60 ( 2000 ).
28. Frontzek , A. V. , Strokov , S. V. , Embs , J. P. & Lushnikov , S. G. Does a dry protein undergo a glass transition? J. Phys. Chem. B 118 , 2796 - 2802 ( 2014 ).
29. Ringe , D. & Petsko , G. A. The 'glass transition' in protein dynamics: What it is, why it occurs, and how to exploit it . Biophys. Chem . 105 , 667 - 680 ( 2003 ).
30. Debenedetti , P. G. & Stillinger , F. H. Supercooled liquids and the glass transition . Nature 410 , 259 - 267 ( 2001 ).
31. Eisenberg , D. , Schwarz , E. , Komaromy , M. & Wall , R. Analysis of membrane and surface protein sequences with the hydrophobic moment plot . J. Mol. Biol . 179 , 125 - 142 ( 1984 ).
32. Bellavia , G. , Giuffrida , S. , Cottone , G. , Cupane , A. & Cordone , L. Protein thermal denaturation and matrix glass transition in different protein-trehalose-water systems . J. Phys. Chem. B 115 , 6340 - 6346 ( 2011 ).
33. Hess , B. , Kutzner , C., van der Spoel , D. & Lindahl , E. Gromacs 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation . J. Chem. Theory Comput . 4 , 435 - 447 ( 2008 ).
34. Van Der Spoel , D. et al. Gromacs: Fast, flexible, and free . J. Comput. Chem . 26 , 1701 - 1718 ( 2005 ).
35. Lindahl , E. , Hess , B. & van der Spoel , D. Gromacs 3.0: A package for molecular simulation and trajectory analysis . J. Mol. Model . 7 , 306 - 317 ( 2001 ).
36. Berendsen , H. J., van der Spoel, D. & van Drunen , R. Gromacs: A message-passing parallel molecular dynamics implementation . Comput. Phys. Commun . 91 , 43 - 56 ( 1995 ).
37. Nos? , S. A unified formulation of the constant temperature molecular dynamics methods . J. Chem. Phys . 81 , 511 - 519 ( 1984 ).
38. Hoover , W. G. Canonical dynamics: Equilibrium phase-space distributions . Phys. Rev. A 31 , 1695 - 1697 ( 1985 ).
39. Parrinello , M. & Rahman , A. Polymorphic transitions in single crystals: A new molecular-dynamics method . J. Appl. Phys . 52 , 7182 - 7190 ( 1981 ).
40. Nos? , S. & Klein , M. Constant pressure molecular dynamics for molecular systems . Mol. Phys . 50 , 1055 - 1076 ( 1983 ).
41. Allen , M. P. & Tildesley , D. J . Computer Simulation of Liquids 64-65 (Oxford University Press, USA, 1989 ).
42. Essmann , U. et al. A smooth particle mesh ewald method . J. Chem. Phys . 103 , 8577 - 8593 ( 1995 ).
43. Hess , B. , Bekker , H. , Berendsen , H. J. C. & Fraaije , J. Lincs: A linear constraint solver for molecular simulations . J. Comput. Chem . 18 , 1463 - 1472 ( 1997 ).
44. Hess , B. P-lincs: A parallel linear constraint solver for molecular simulation . J. Chem. Theory Comput . 4 , 116 - 122 ( 2008 ).
45. Miyamoto , S. & Kollman , P. A. Settle : An analytical version of the shake and rattle algorithm for rigid water models . J. Comput. Chem . 13 , 952 - 962 ( 1992 ).
46. Best , R. B. & Mittal , J. Protein simulations with an optimized water model: Cooperative helix formation and temperature-induced unfolded state collapse . J. Phys. Chem. B 114 , 14916 - 14923 ( 2010 ).
47. Duan , Y. et al. A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations . J. Comput. Chem . 24 , 1999 - 2012 ( 2003 ).
48. Abascal , J. L. F. & Vega , C. A general purpose model for the condensed phases of water: Tip4p/2005 . J. Chem . Phys. 123 , 234505 ( 2005 ).
49. Flyvbjerg , H. & Petersen , H. G. Error estimates on averages of correlated data . J. Chem. Phys . 91 , 461 - 466 ( 1989 ).
50. Grossfield , A. & Zuckerman , D. M. Quantifying uncertainty and sampling quality in biomolecular simulations . Annu. Rep. Comput. Chem . 5 , 23 - 48 ( 2009 ).