Microcanonical thermostatistics analysis without histograms: cumulative distribution and Bayesian approaches
Microcanonical thermostatistics analysis has become an important tool to reveal essential aspects of phase transitions in complex systems. An efficient way to estimate the microcanonical inverse temperature and the microcanonical entropy is achieved with the statistical temperature weighted histogram analysis method (ST-WHAM). The strength of this method lies on its flexibility, as it can be used to analyse data produced by algorithms with generalised sampling weights. However, for any sampling weight, ST-WHAM requires the calculation of derivatives of energy histograms , which leads to non-trivial and tedious binning tasks for models with continuous energy spectrum such as those for biomolecular and colloidal systems. Here, we discuss two alternative methods that avoid the need for such energy binning to obtain continuous estimates for in order to evaluate by using ST-WHAM: (i) a series expansion to estimate probability densities from the empirical cumulative distribution function (CDF), and (ii) a Bayesian approach to model this CDF. Comparison with a simple linear regression method is also carried out. The performance of these approaches is evaluated considering coarse-grained protein models for folding and peptide aggregation.
pacs:07.05.Kf, 02.50.Ng, 02.60.Jh
Fundamental aspects of phase transitions in complex systems can be revealed by the analysis of its microcanonical thermostatistics gross-book (); gross2 (), which is characterised by the well known entropy , where denotes the density of states of a system with energy . In particular, the analysis of inflection points of the microcanonical inverse temperature plays an important role in the identification of stable, unstable and metastable regions in the phase diagram frigori2010JCS (); schnabelPRE2011 (); rochaPRE2014 (), providing alternative insights to the usual canonical analysis. Also, free-energy profiles can be obtained from the caloric curves vs , from where one can easily evaluate the values of barrier heights and latent heats. In this way, the microcanonical thermostatistics analysis has been incorporated in many studies in the literature, e.g. Refs. frigoriJCP2013 (); jankePRL2006 (); jankeJCP2008 (); moddel2010pccp (); bereau2010jacs (); church2012jcp (); straubPRL2014 () to name a few.
Nevertheless, any analysis must rely on data obtained from efficient exploration of the configurational space. It is well known that numerical simulations performed with Monte Carlo (MC) and molecular dynamics (MD) methods pose limitations to the achievement of reliable data sampling okamoto-60 (). Such limitations are related to the critical slowing down landau-book (), which is observed in studies of continuous phase transitions, and to the entrapment in local minima, in the case of systems with rugged energy landscapes. In both cases, the configurational space is poorly explored in a reasonable computational simulation time, which may produce biased physical averages. To overcome the trapping problem, it has been suggested that configurations must be sampled using algorithms based on generalised ensembles, where the updates are performed with non-Boltzmann statistical weights . For instance, the multicanonical algorithm (MUCA) berg2000 (); berg2003 (), the extended Gaussian ensemble (EGE) johal (); neuhaus06 (); frigoriEurP2010 (), Tsallis statistical weight straub-tsallis (); hansmannPRE1997 (), and replica exchange method (REM) hukushima1996 () either use a series of Boltzmann weights or any convenient generalised sampling weight straub-gREM ().
MUCA simulations sample configurations with a weight in such a way that the energy distribution is uniform, . Thus, a precise determination of is equivalent to obtain an estimate for the density of states , i.e. . The weights follows from the parameterization of the entropy , where and are the so-called multicanonical parameters. The iterative procedure to obtain the MUCA parameters is described in detail in references berg2000 (); berg2003 (), and can be read as,
where stands for the th multicanonical simulation. The recursion steps require the discretisation of the energy for continuous energy models. Therefore, it is convenient to define , where is the binsize, is an integer, and is a constant that defines a reference energy. All the energies in the interval contribute to the histogram .
Methods to improve sampling based on simulations at different temperatures have been proposed to either be conducted in parallel (REM) or as a random walk between different temperatures. In REM, non-interacting replicas of the system are simultaneously simulated by the usual MC or MD algorithms, and from time to time, pairs of replicas at neighboring temperatures are exchanged with a transition probability. From the data produced by simulations performed at a single temperature or at a set of temperatures , with , it is necessary to employ a reweighing scheme to evaluate physical averages at a given temperature . Reweighting techniques ferrenberg-rewei (); wham (); alves92 () use data from either a single histogram or multiple histograms obtained from MC or MD simulations.
Recently, a simple method called statistical weighted histogram analysis method (ST-WHAM) kimJCP-stwham () has been proposed as an iteration-free procedure to obtain an estimate for the microcanonical inverse temperature. In this method the usual WHAM equations ferrenberg-rewei (); wham () are converted into a weighted average of the individual densities of states obtained from simulations carried out with different sampling weights . From energy histograms produced by multiple simulations, ST-WHAM yields a statistical temperature , which is an estimate of the inverse microcanonical temperature . Interestingly, there is a numerical procedure based on the multicanonical recursion relations (1) and (2), which is called ST-WHAM-MUCA rizziJCP2011 (), that can be used to replace the direct integration in order to evaluate the entropy . Although both ST-WHAM and ST-WHAM-MUCA have the advantage of a posteriori discretisation of energies, their naive implementations may lead to biased evaluations of physical quantities for continuous energy models just like all the aforementioned algorithms.
As described in Section II, the estimates for inverse microcanonical temperature depends on the derivatives of the energy histograms (see Eq. (4)). Here, we analyse how the estimates are energy binning dependent and, in Section III, we present two alternative approaches that avoids the need for energy binning to evaluate the microcanonical caloric curve for continuous energy models: (i) a proposal by Berg and Harris berg-accum (), which involves an empirical cumulative distribution (CDF) and uses discrete Fourier series; and (ii) a Bayesian approach bayes-caldwell () to model this CDF and assumes that the thermodynamic phase transitions are well described by the coexistence of two phases. A comparative analysis between these approaches is made in order to characterise for two systems that undergo first-order-like phase transitions: the folding transition of a coarse-grained protein model for Ubiquitin and the aggregation transition of two heteropolymers that follows a Fibonacci sequence. These examples allow us to describe the statistical and systematic errors involved in the numerical calculations of and , which are presented in Section IV. Conclusions on this comparative analysis are presented in Section V.
Ii Statistical temperature weighted histogram method
The ST-WHAM kimJCP-stwham () yields a direct estimate of the inverse microcanonical temperature by considering the statistical inverse temperature
where , and . It is preferable to write Eq. (3) as
This ST-WHAM-MUCA algorithm is quite simple if one has .
Iii Numerical Evaluation of Derivatives
iii.1 Linear regression
We can numerically evaluate the derivatives in Eq. (4) in a naive way, where the derivatives at energies follow from a linear regression around this point. For instance, we use a linear regression with points; selecting points means that the derivative at is calculated with the values of , where . We chose a value for according to the energy binsize . Consequently, we calculate the derivatives in the energy range . In this method, it is more convenient to directly calculate the derivative of than the derivative of . We calculate the linear regression with a subroutine easily adapted from the linear fit subroutine in recipes ().
iii.2 Cumulative distribution method
Another approach can be devised by considering an algorithm based on the cumulative distribution function (CDF) berg-accum (). The advantage of such approach is that it avoids histogramming when describing probability densities , dismissing the need for any ad hoc energy discretisation. The method defines an estimator for the CDF , where the function is an empirical cumulative distribution function (ECDF) for the probability density . The algorithm sorts the energy time series of length NDAT in an ascending order , so any outliers can be eliminated by constructing a restricted ECDF in the range between two meaningful points and (in general one takes and ). The basic idea is to propose an approximating function to describe , from where the difference function is defined,
This function can be expanded in Fourier series,
which gives the Fourier coefficients berg-accum (),
Notably, Eq. (8) provides all the information that one needs to obtain a parameter-free energy probability density . The reason is because Eq. (6) yields a smooth estimate of the CDF given MMAX coefficients in the Fourier expansion,
even if one assumes a linear ansatz for (see berg-accum ()). In turn, it becomes easy to obtain estimates for the probability density by differentiation, which consists in a smooth estimation of .
Now, let us go back to the numerical differentiation in Eq. (4). To obtain parameter-free derivatives , we just need to differentiate again. This corresponds to the following numerical calculation,
Unlike the linear regression method, where we calculate the derivative of , here is more convenient to compute the derivative of directly. We have included in Appendix the FUNCTION DERPD(X), which can be used to estimate this derivative from the Fourier coefficients . One can easily incorporate this function in the program CDF_PD subroutine provided by Berg and Harris berg-accum ().
iii.3 Bayesian analysis
Next we present our approach to implement a Bayesian analysis. To this end, input data (denoted by ) is the empirical cumulative distribution and not the histograms that are dependent on the energy binsize . Thus, the aim is to describe the empirical cumulative distribution by a given model (called model M) and using the Bayesian analysis bayes-caldwell (), extract values for its parameters (denoted by ) from the posterior probability density function (PDF) . Here, denotes possible parameters related to the experimental conditions, which is the temperature in our numerical experiment, and refers to a set of observed experimental points . From the posterior PDF, we extract the desired probability distributions for each parameter,
which correspond to the marginalized distributions. These marginalized distributions are depicted with histograms in Fig. 6 for the Ubiquitin.
For each protein, our experimental data is produced by dividing the observed energy range into 35 points and experimental errors are extracted from the jackknife procedure. We anticipate here that these experimental data are displayed in Figures 5 and 11 for the folding-unfolding transition of Ubiquitin and for the aggregation of two Fibonacci sequences, respectively. Within a Bayesian approach, the set of observations is just a set among possible experimental outcomes at the points . Thus, given a fixed model M and a set of experimental values , parameters estimation within Bayesian approach is obtained as follows bayes-caldwell (),
where is the prior distribution for a fixed model M.
Since our aim is to characterise in a first-order phase transition region, we assume that the energy distribution (i.e. a continuous estimate for in Eq. (4)) is well described by a double-peak function. This function characterises the coexistent phases at temperatures close enough to the transition temperature , whose shape is a consequence of the free-energy barrier between the ordered and disordered phases. Therefore, the expected energy distribution can be written as a normalized sum of two Gaussian distributions with peaks centered at energies and , corresponding to the disordered and ordered phases, respectively Lee-1991 (). If we assume the model M as the cumulative distribution function of , and recalling that for a Gaussian distribution the CDF is the error function, we arrive to the following modelling,
This is a 5-parameter model, which we call model I in our analysis, with .
Now, to obtain the posterior probability of the parameter set , given data set , from the likelihood times the initial prior probabilities for our model I, we use a Markov chain Monte Carlo. To this end, we consider the Metropolis algorithm with the proposal function
where are the errors in our computational experiment. The values are calculated with the accepted values for the parameters of model I.
Iv Numerical Comparisons
In the following, we present results where we evaluate the performance of the different approaches in obtaining continuous estimates for and, consequently, the estimate for microcanonical inverse temperature. The comparative analysis is made with data obtained from both MC and MD simulations, where we employ two simplified off-lattice protein models, i.e. all atoms in the polypeptide chain are replaced by beads located at the -atom position, and present a continuous energy spectrum.
Ubiquitin is a 76 residue small protein (PDB code 1UBQ) for which evidences suggest a two-state folding mechanism. It had received considerable attention for what concerns solvent response and temperature dependence in the secondary structure formation ubqJacob (); chung-ubiquitin (); ghosh-ubiquitin () and stretching experiments jernigan-ubiquitin (); das-ubiquitin (); imparato-ubiquitin ().
To study the performance of the numerical approaches in evaluating through the folding-unfolding transition, we perform simulations with a coarse-grained version for Ubiquitin but with a rather detailed potential. We use the structure-based model implemented in SMOG@ctbp web server smog () to perform MD simulations with GROMACS package. The analyses were made from measurements during 250 ns of simulation after a period of thermalization. Measurements were obtained from 7 independent simulations at temperatures with and .
First, we evaluate the performance of the approaches considering the data obtained from a single MD trajectory produced at , which is close to the transition temperature evaluated with statistics obtained with 7 independent MD simulations. Figure 2 displays the estimates for the caloric curve using the linear regression to evaluate derivatives of from data obtained at . This figure shows how noisy the caloric curve can be as a function of the energy binsize and 1.3, keeping the number of points fixed. Figure 2 compares what seems a good energy discretisation with the CDF method. Error bars based on 20 jackknife bins are presented only for the CDF method. Now, Fig. 3 compares both methods when we consider the statistics obtained with the entire set of temperatures , . To present error bars for the CDF method, we have increased the number of jackknife bins to 80 because this larger statistics. Nice agreement between both methods are obtained if we choose conveniently.
Figure 5 displays what we call input experimental points to perform the Bayesian analysis. This experimental data was obtained from MD simulations at . The dashed line corresponds to the cumulative distribution obtained with model I. The parameters of this model that fit the input data are displayed in Table 1. We include values calculated from means of the marginalized distributions, and global modes that presented the smallest /d.o.f of the model. To verify how the proposed model describes the energy distribution when compared with the CDF method and naive histogramming, we show in Fig. 5 the results obtained from these methods for the statistics collected at . We realize that model I does not recover properly the region between the peaks of the energy distributions obtained with previous numerical methods. The inset in this figure compares the calculations of following from the CDF and Bayesian approaches. In Fig. 6 we display the marginalized distributions of , including the correlation between the parameters and to illustrate their interdependence. We have always used flat priors over appropriated ranges to obtain the parameter distributions.
iv.2 Fibonacci sequences
Figure 8 shows the caloric curve for the aggregation of two chains of monomers described by the AB model, a coarse-grained off-lattice protein model that replaces the all-atom potentials by simpler physical interactions. This model reduces the protein to a chain of monomers of two types: the hydrophobic () and the polar or hydrophilic () types, located at the atoms. The AB model has also been used in studies to explore aggregation phenomenon jankePRL2006 (); jankeJCP2008 (), where more than one chain is included in the system, and in studies of microcanonical thermostatistics of heteropolymers frigoriJCP2013 ().
Here, we consider a simple system which consists of two heteropolymers defined by Fibonacci sequences with monomers, i.e. ABBABBABABBAB. The statistics for this system amounts to sweeps per replica produced with REM. Attempts to exchange the replicas occur every sweeps, using a scheme that alternates attempts between even replicas and their neighbors, and odd replicas and their neighbors. Although one has several manners to define the temperature set fiore (), here we determine it using an arithmetic progression . We consider replicas choosing , and .
Figure 8 depicts the caloric curves obtained with energy binsizes , and . The value was used to calculate the derivatives of . Different values of allows one to verify the systematic errors associated with and the use of linear regression. For a fixed , small values of (or ) tend to reproduce the statistical fluctuations of the data. On the other hand, a value as large as eliminates the physical information associated with the small S-loop at . For comparison, we include in Fig. 8 the estimates from the CDF method. These comparisons are based on a single time series selected from data produced with REM at the inverse temperature . It is important to reliably calculate the caloric curve when assessing the canonical critical temperature and the latent heat of the transition. We can figure out how important the fluctuations observed in these curves are by calculating the statistical errors associated with the CDF method. Figure 8 displays the statistical error bars calculated with 20 jackknife bins for the dataset with measurements obtained with REM at . Both methods present comparable statistical error bars. This figure shows that the small S-loop at does not result from the statistical fluctuations in our dataset. Therefore, any smooth estimate of this curve at , obtained for example with (), would hide physical information.
In Fig. 9 we include data produced with replicas together. Comparison of the results presented in this figure with those obtained by multicanonical simulations in Refs. jankePRL2006 (); jankeJCP2008 () reveals that they are in excellent agreement: even the small S-loop near appears in our case. Computation of the inverse of aggregation temperature also resulted in a quantitative agreement. Larger error bars appear around , probably because the lack of temperatures in the set used to exchange the replicas. The averaged estimates from both derivative calculation methods agreed for a particular choice of and points in the linear regression method.
Now, Fig. 11 shows the experimental data we have used to perform the Bayesian analysis. Again, we initially keep the analysis restricted to a single energy time series, in this case produced at . We have included the Bayesian estimates (dashed line) of the cumulative function assuming the model I. The parameters for this model are presented in Table I. To verify how well this model reproduces , we include in Fig. 11 the energy distributions obtained with a naive histogramming procedure, CDF method and Bayesian analysis. Clearly, Bayesian analysis does not reproduce the expected behaviour between both peaks. This leads to a biased estimation of as can be seem in the inset of this figure. Moreover, the approach with model I does not reproduce the small S-loop observed at (data not shown), in contrast to the estimates with a convenient choice of or the CDF method. As a consequence, we again do not follow further analysis with the Bayesian approach including all statistics from the 7 replicas.
V Summary and conclusions
We presented two alternative approaches for the estimation of the energy probability distributions avoiding the need of energy binning in order to obtain the microcanonical thermostatistics analysis from ST-WHAM. Numerical comparisons between the approaches were presented and we showed the statistical and systematic errors that can arise when one evaluates the microcanonical inverse temperature for two continuous energy models that exhibits first-order like phase transitions. Our results indicate that the CDF method yields reliable estimates for both and when compared with the linear regression method, and far more successfully than the Bayesian approach with model I. Unlike the linear regression method, the CDF method avoids the undesirable task of choosing the energy binsize for a careful evaluation of the microcanonical temperature, as highlighted in the analysis of the aggregation transition of the two Fibonacci sequences. In this case, the caloric curve presents a quite unusual behaviour with two S-loops, indicating two transitions. In particular, we showed that the use of large values for in linear regression method would hide physical information about the small S-loop at . On the other hand, the small S-loop could not be detected with the Bayesian analysis (data not shown) as a consequence of the poor evaluation of and its derivative (Fig. 11). The reason is because the model I consists of only 5 parameters. This is a low number compared to the CDF method, which allows a variable number of parameters (i.e. Fourier coefficients) defined depending on the demand of the ECDF. For instance, the CDF method needs 13 Fourier coefficients to obtain probability densities for the Ubiquitin data and this number goes up to 74 for Fibonacci sequences. Furthermore, model I was constructed on the hypothesis that the coexistence of two thermodynamic bulk phases can be well described by two gaussian distributions. Actually, this is an approximation because the energies in between the two peaks, describing mixed phase configurations, are not properly take into account in the two gaussian model.
We thank Bernd Berg for carefully reading our manuscript. This work has been supported by the Brazilian agencies FAPESP, CNPq and RUSP (University of São Paulo).
Derivative of the probability density:
FUNCTION DERPD(X) include ’../../Libs/Fortran/implicit.sta’ include ’../../Libs/Fortran/constants.par’ PARAMETER(NMAX=100) COMMON /CDFProb/ XMIN,XRANGE,DN(NMAX),M !M number of Fourier coefficients DERPD=ZERO DO J=1,M AUX=J*PI/XRANGE AUX=AUX*AUX DERPD=DERPD-DN(J)*AUX*SIN(J*PI/XRANGE*(X-XMIN)) ENDDO RETURN END
- (1) D.H.E. Gross, Microcanonical thermodynamics: phase transitions in small systems, Lecture Notes in Physics 66 (World Scientific, Singapore, 2001).
- (2) D.H. Gross, and J.F. Kenney, J. Chem. Phys. 122 (2005) 224111.
- (3) R.B. Frigori, L.G. Rizzi, and N.A. Alves, J. Phys.: Conf. Ser. 246 (2010) 012018.
- (4) S. Schnabel, D.T. Seaton, D.P. Landau, and M. Bachmann, Phys. Rev. E 84 (2011) 011127.
- (5) J.C.S. Rocha, S. Schnabel, D.P. Landau, and M. Bachmann, Phys. Rev. E 90 (2014) 022601.
- (6) R.B. Frigori, L.G. Rizzi, and N.A. Alves, J. Chem. Phys. 138 (2013) 015102.
- (7) C. Junghans, M. Bachmann, and W. Janke, Phys. Rev. Lett. 97 (2006) 218103.
- (8) C. Junghans, M. Bachmann, H. Arkin, and W. Janke, J. Chem Phys. 128 (2008) 085103.
- (9) M. Möddel, W. Janke, and M. Bachmann, Phys. Chem. Chem. Phys. 12 (2010) 11548.
- (10) T. Bereau, M. Bachmann, and M. Deserno, J. Am. Chem. Soc. 132 (2010) 13129.
- (11) M.S. Church, C.E. Ferry and A. E. van Giessen, J. Chem. Phys. 136 (2012) 245102.
- (12) W.J. Cho, J. Kim, J. Lee, T. Keyes, J.E. Straub, and K.S. Kim, Phys. Rev. Lett. 112 (2014) 157802.
- (13) A. Mitsutake, Y. Sugita, and Y. Okamoto, Biopolymers 60 (2001) 96.
- (14) D.P. Landau, and K. Binder, A guide to Monte Carlo simulations in Statistics Physics. Cambridge University Press, 2000.
- (15) B.A. Berg, 2000. Monte Carlo Methods, Fields Institute Communications Vol. 26, edited by N. Madras (American Mathematical Society, Providence, RI, 2000), p. 1-24.
- (16) B.A. Berg, Comput. Phys. Commun. 153 (2003) 397.
- (17) R.S. Johal, and A. Planes, E. Vives, Phys. Rev. E 68 (2003) 056113.
- (18) T. Neuhaus, and J.S. Hager, Phys. Rev. E 74 (2006) 036702.
- (19) R.B. Frigori, L.G. Rizzi, and N.A. Alves, Eur. Phys. J. B 75 (2010) 311.
- (20) I. Andricioaei, and J.E. Straub, Phys. Rev. E 53 (1996) R3055.
- (21) U.H.E. Hansmann, and Y. Okamoto, Phys. Rev. E 56 (1997) 2228.
- (22) K. Hukushima, and K. Nemoto, J. Phys. Soc. Japan 65 (1996) 1604.
- (23) J. Kim, T. Keyes, and J.E. Straub, J. Chem. Phys. 132 (2010) 224107.
- (24) A.M. Ferrenberg, and R.H. Swendsen, Phys. Rev. Lett. 61 (1988) 2635. [Erratum: 63 (1989) 1658]
- (25) A.M. Ferrenberg, and R.H. Swendsen, Phys. Rev. Lett. 63 (1989) 1195.
- (26) N.A. Alves, B.A. Berg, and S. Sanielevici, Nucl. Phys. B 376 (1992) 218.
- (27) J. Kim, T. Keyes, and J.E. Straub, J. Chem. Phys. 135 (2011) 061103.
- (28) L.G. Rizzi, and N.A. Alves, J. Chem. Phys. 135 (2011) 141101.
- (29) B.A. Berg, and R.C. Harris, Comput. Phys. Commun. 179 (2008) 443.
- (30) A. Caldwell, D. Kollár, and K. Kröninger, Comput. Phys. Comm. 180 (2009) 2197.
- (31) W. Press et al., Numerical Recipes (Cambridge University Press, London, 1986).
- (32) J. Lee, and J.M. Kosterlitz, Phys. Rev. B 43 (1991) 3265.
- (33) J. Jacob, B. Krantz, R.S. Dothager, P. Thiyagarajan, and T.R. Sosnick, J. Mol. Biol. 338 (2004) 369.
- (34) H.S. Chung, and A. Tokmakoff, Proteins 72 (2008) 474.
- (35) S.G. Dastidar, and C. Mukhopadhyay, Phys. Rev. E 72 (2005) 051928.
- (36) J.I. Sulkowska, A. Kloczkowski, T.Z. Sen, M. Cieplak, R.L. Jernigan, Proteins 71 (2008) 45.
- (37) A. Das, C. Mukhopadhyay, Proteins 75 (2009) 1024.
- (38) A. Imparato, A. Pelizzola, Phys. Rev. Lett. 100 (2008) 158104.
- (39) J.K. Noel, P.C. Whitford, K.Y. Sanbonmatsu, and J.N. Onuchic, Nucl. Acids Res. 38 (2010) W657.
- (40) C.E. Fiore, J. Chem. Phys. 135 (2011) 114107.