Reweighting distributions between thermodynamic states

Frank Noefrank.noe@fu-berlin.de
FU Berlin, Arnimallee 6, 14195 Berlin, Germany

Summary: Reweighting is a principle used to analyze simulations at different thermodynamic states, such as temperatures or Hamiltonians (e.g. parallel tempering or umbrella sampling). The reweighting principle states that at each thermodynamic state a configuration $\mathbf{x}$ has certain equilibrium distribution, and the equilibrium distributions at different thermodynamic states can be reweighted into on another by multiplying the density at one thermodynamic state with a reweighting factor that depends only on $\mathbf{x}$ and the two states to be related. This principle is used by various estimators, including the weighted histogram analysis method (WHAM), Bennett’s acceptance ratio (BAR), the multistate Bennett acceptance ratio (MBAR) and the transition-based reweighting analysis method (TRAM).


Enhanced sampling simulations to which the reweighting principle applies (at least partially) include:

  • Unbiased molecular dynamics or Metropolis (Hastings) Monte Carlo simulations

  • Generalized ensemble methods, such as replica exchange molecular dynamics [6], parallel tempering [10][9][21] and simulated tempering [14]

  • Umbrella sampling [23] which uses a set of biased Hamiltonians to ensure approximately uniform sampling along a set of pre-defined slow coordinates.

  • Metadynamics and conformational flooding [13][7].

  • Accelerated dynamics [8]

Reduced potential and statistical weight

A configuration $\mathbf{x}$ is a point in configuration space $\Omega$ containing all quantities characterizing the instantaneous state of the system, such as particle positions or spins, system volume (in constant-pressure ensembles), and particle numbers (in ensembles of constant chemical potential).

We consider a set of simulation trajectories, $\{\mathbf{x}^{I}(t)\}$, sampling from $\Omega$, each at a given thermodynamic state $I$. A thermodynamic state, here denoted as captial superscript letters $I$, $J$, $K$, is characterized by its thermodynamic variables, such as temperature, pressure or chemical potential. The dynamics are assumed to fulfill microscopic detailed balance at their respective thermodynamic states.

We follow the notation of [19] and write a generalized potential function. In a thermodynamic state $I$, defined by a particular combination of the potential energy function $U^{I}$, pressure $p^{I}$, chemical potentials $\mu_{s}^{I}$ of chemical species $s$ and temperature $T^{I}$, our system has a reduced potential defined by: $$u^{I}(\mathbf{x})=\frac{U^{I}(\mathbf{x})+p^{I}V(\mathbf{x})+\sum_{s}\mu_{s}^{I}N_{s}(\mathbf{x})}{k_{B}T^{I}}. \:\:\:\:(1)$$ Here, $V(\mathbf{x})$ is the volume of the system in configuration $\mathbf{x}$ and $N_{s}(\mathbf{x})$ count the particle numbers of species $s$ at configuration $\mathbf{x}$.

The probability density of configuration $\mathbf{x}$ can, for any arbitrarily chosen thermodynamic state, be expressed as: $$\rho^{I}(\mathbf{x})=\frac{1}{Z^{I}}e^{-u^{I}(\mathbf{x})}, \:\:\:\:(2)$$ where $Z^{I}$ is the partition function of thermodynamic state $I$:

$$Z^{I}=\int_{\Omega}d\mathbf{x}\: e^{-u^{I}(\mathbf{x})}. \:\:\:\:(3)$$

Reweighting between thermodynamic states

The principle of reweighting is the following: consider a reference thermodynamic state $I$, and a different thermodynamic state $J$. We want to express the equilibrium probability in $J$ with respect to the equilibrium probability in $I$. This can be expressed as: $$\begin{aligned} \rho^{J}(\mathbf{x}) & = & \frac{1}{Z^{J}}e^{-u^{J}(\mathbf{x})}\nonumber \\ & = & \frac{1}{Z^{J}}e^{-(u^{J}(\mathbf{x})-u^{I}(\mathbf{x}))}e^{-u^{I}(\mathbf{x})}\nonumber \\ & = & \frac{1}{Z^{J}}e^{-b^{J}(\mathbf{x})}e^{-u^{I}(\mathbf{x})} \:\:\:\:(4)\end{aligned}$$ and equivalently as $$\rho^{J}(\mathbf{x})=\frac{Z^{I}}{Z^{J}}\gamma^{J}(\mathbf{x})\rho^{I}(\mathbf{x}) \:\:\:\:(5)$$ Here, the energy difference $$b^{J}(\mathbf{x})=u^{J}(\mathbf{x})-u^{I}(\mathbf{x}) \:\:\:\:(6)$$ is the bias energy of thermodynamic state $J$ with respect to the reference state $I$, and $$\gamma^{J}(\mathbf{x})=e^{-b^{J}(\mathbf{x})} \:\:\:\:(7)$$ is the corresponding reweighting factor. Eq. (5) shows that density at state $J$ is proportional to the density at state $I$, multiplied by a reweighting factor that can be expressed using a bias energy $b^{J}(\mathbf{x})$. Equations (4-5) will be called point-wise reweighting, as they reweigh individual sample configurations and do not require a discretization of reweighting energies or factors.

The reweighting formulation is useful because in the aforementioned enhanced sampling simulations the bias energy $b^{J}(\mathbf{x})$ is known. We can thus simulate in biased ensembles, and by virtue of Eq. (5) reweight our sampled distribution $\rho^{J}(\mathbf{x})$ to the distribution of interest, $\rho^{I}(\mathbf{x})$.

Example 1: umbrella sampling. When the simulation setup contains multiple biased simulations, such as in umbrella sampling or metadynamics, there is usually a unique temperature $T$, but different potentials $U^{I}(\mathbf{x})=U(\mathbf{x})+B^{I}(\mathbf{x})$ are employed where $B^{I}(\mathbf{x})$ is the $I$th bias potential. The original potential is given by $$u(\mathbf{x})=\frac{U(\mathbf{x})}{k_{B}T}. \:\:\:\:(8)$$ with a corresponding statistical weight $\rho(\mathbf{x})\propto\mathrm{e}^{-u(\mathbf{x})}$. The biased potential is given by: $$u^{I}(\mathbf{x})=\frac{U(\mathbf{x})+B^{I}(\mathbf{x})}{k_{B}T}. \:\:\:\:(9)$$ We can reweigh the unbiased to the biased distributions via Eq. (4), using $$b^{I}(\mathbf{x})=\frac{B^{I}(\mathbf{x})}{k_{B}T} \:\:\:\:(10)$$ as bias energy.

Example 2: parallel tempering. In PT or REMD simulations, the weighting occurs through the different temperatures: Given a configuration $\mathbf{x}$ with potential energy $U(\mathbf{x})$, the statistical weight at temperature $T$ is proportional to $\mathrm{e}^{-u(\mathbf{x})}$ using the reduced potential energy $$u(\mathbf{x})=\frac{U(\mathbf{x})}{k_{B}T} \:\:\:\:(11)$$ at another temperature $T^{I}$, we have $$u^{I}(\mathbf{x})=\frac{U(\mathbf{x})}{k_{B}T^{I}} \:\:\:\:(12)$$ and statustical weight $\rho^{I}(\mathbf{x})\propto\mathrm{e}^{-u^{I}(\mathbf{x})}$. Reweighting via Eq. (4) occurs as: $$\begin{aligned} \rho^{I}(\mathbf{x}) & = & \frac{1}{Z^{J}}e^{-(u^{I}(\mathbf{x})-u(\mathbf{x}))}e^{-u(\mathbf{x})}\\ & = & \frac{1}{Z^{J}}e^{-\frac{U(\mathbf{x})}{k_{B}T^{*}}}e^{-u(\mathbf{x})}\\ T^{*} & = & \frac{T\, T^{I}}{T^{I}-T} \:\:\:\:(13)\end{aligned}$$ thus the reweighting factor and also the bias energy $$b^{I}(\mathbf{x})=\frac{U(\mathbf{x})}{k_{B}T^{*}} \:\:\:\:(14)$$ depend on the potential energy of the system.

Discrete reweighting

In special cases the state space $\Omega$ is discrete (e.g. when our dynamics is given by a Markov chain or Markov jump process). In some other cases, $\Omega$ is continuous, but the state space can be discretized such that the reweighting factors $\gamma^{J}(\mathbf{x})$ are well represented by a discrete approximation: $$\gamma^{J}(\mathbf{x})\approx\gamma_{i}^{J}=\gamma^{J}(\mathbf{x}_{i}) \:\:\:\:(15)$$ where $\mathbf{x}_{i}$ is a representative point close to $\mathbf{x}$. An example is umbrella sampling, where the bias energy is usually a function of only a few (often one) order parameters $f_{1}(\mathbf{x}),\, f_{2}(\mathbf{x}),...$. In this case we can efficiently discretize these order parameters with a fine grid and therefore obtain a good discretization of $\gamma^{J}(\mathbf{x})$ despite the exponential dependence on the bias energy $b^{J}(\mathbf{x})$.

If such a discretization can be made, we call the unbiased, or reference probability, and the biased probability at thermodynamic state $k$: $$\begin{aligned} \pi_{i} & = & \int_{s_{i}}\rho(\mathbf{x})\: d\mathbf{x}\\ \pi_{i}^{(k)} & = & \int_{s_{i}}\rho^{(k)}(\mathbf{x})\: d\mathbf{x} \:\:\:\:(16)\end{aligned}$$ Analogous to point-wise reweighting, they are related by a known and constant bias factor $$\gamma_{i}^{(k)}=\mathrm{e}^{-b_{i}^{(k)}} \:\:\:\:(17)$$ yielding the discrete reweighting equation: $$\pi_{i}^{(k)}=f^{(k)}\pi_{i}\gamma_{i}^{(k)}, \:\:\:\:(18)$$ with the normalization constant $$f^{(k)}=\frac{1}{\sum_{l}\pi_{l}\gamma_{l}^{(k)}} \:\:\:\:(19)$$ Thus, the bias is multiplicative in probabilities or additive in the potential. Note that this discrete reweighting formulation assumes that $b_{i}^{(k)}$ is constant over set $S_{i}$, which requires a very fine discretization.

Unfortunately, this discretization approach does not work well for temperature-based methods because for such method - see Eq. (14) - the reweighting factor depends on the potential energy of the system, which is a property that has a very large range for many-body systems. The standard deviation of the potential energy distribution is proportional to $\sqrt{d}$, where $d$ is the number of degrees of freedom. In typical simulation systems, the distribution of $U(\mathbf{x})$ can cover 100s to 1000s of $kT$. At least 10 bins per $kT$ are needed in order to arrive at a discretization that is acceptable for discrete reweighting. Since the potential energy is usually not the quantity of interest to the analyst, this discretization must be combined with yet another discretization in other degrees of freedom. These requirements typically lead to an inacceptable number of discretization bins that would conterminate the analysis with statistical errors. Thus, with temperature-based sampling methods, point-wise reweighting should be applied, using Eq. (5).

Reweighting estimators

Usually we will have simulations at not only two but at thermodynamic states, and we want to obtain an optimal estimate for the equilibrium distribution at a state of intereste using all the data. This is exactly the purpose of reweighting estimators. The most prominent examples are:

  1. The weighted histogram analysis method (WHAM) [5][12] performs an optimal reweighting amongst multiple thermodynamic states for discrete reweighting factors.

  2. Bennett’s acceptance ratio (BAR) [2] performs optimal point-wise reweighting by virtue of Eq. (5) between two thermodynamic states.

  3. Multistate Bennett acceptance ratio (MBAR) [1][11][19], also known as binless WHAM [20][22] generalizes BAR to multiple thermodynamic states. In the case of discrete state spaces or when the reweighting factors can be efficiently discretized, MBAR is identical to WHAM

  4. The discrete transition-based reweighting analysis method (dTRAM) [25][24] is a generalization to both WHAM and reversible Markov models. It is a discrete reweighting estimator that also takes the time sequence of trajectories into account. dTRAM can obtain superior estimates of the equilibrium distribution compared to WHAM because taking the time correlation in the data into account allows for an increased statistical efficiency. Additionally, dTRAM can provide estimates of the kinetics at all thermodynamic states.

  5. Dynamic histogram analysis method (DHAM) [18] addresses a similar problem like dTRAM and employs diffusional dynamics as a model to facilitate a solution of the likelihood problem.

  6. The continuous transition-based reweighting analysis method (TRAM) [15] is a generalization to both MBAR and reversible Markov models. It is a point-wise reweighting estimator that also takes the time sequence of trajectories into account. TRAM can obtain superior estimates of the equilibrium distribution compared to MBAR because taking the time correlation in the data into account allows for an increased statistical efficiency. Additionally, TRAM can provide estimates of the kinetics at all thermodynamic states.

  7. Dynamical reweighting [16][3][17][4] applies the idea of reweighting to the probabilities of trajectory pieces rather than configurations, and can therefore employ generalized ensemble or biased simulations in order to estimate dynamical expectation values.

Acknowledgements:

Parts of this article have been published in [24] and [15].

Bibliography

[1]: C. Bartels: Analyzing biased Monte Carlo and molecular dynamics simulations. Chem. Phys Lett. 331, 446-454 (2000).

[2]: Charles H. Bennett: Efficient Estimation of Free Energy Differences from Monte Carlo Data. J. Chem. Phys. 22, 245-268 (1976).

[3]: John D. Chodera, William C. Swope, Frank Noé, Jan-Hendrik Prinz and Vijay S. Pande: Dynamical reweighting: Improved estimates of dynamical properties from simulations at multiple temperatures. J. Phys. Chem. 134, 244107 (2011).

[4]: Du, Wei-Na, Marino, Kristen A. and Bolhuis, Peter G.: Multiple state transition interface sampling of alanine dipeptide in explicit solvent. J. Chem. Phys. 135, 145102 (2011).

[5]: Alan M. Ferrenberg and Robert H. Swendsen: Optimized Monte Carlo data analysis. Phys. Rev. Lett. 63, 1195-1198 (1989).

[6]: C. G. Greyer: Computig science and statistics, the 23rd symposium of the interface (Interface Foundation, Fairfax). American Statistical Association, New York. (1991).

[7]: Grubmüller, H.: Predicting slow structural transitions in macromolecular systems: conformational flooding. Phys. Rev. E 52, 2893 (1995).

[8]: D. Hamelberg, J. Mongan and J. A. McCammon: Accelerated molecular dynamics: A promising and efficient simulation method for biomolecules. J. Chem. Phys. 120, 11919-11929 (2004).

[9]: Hansmann, Ulrich H.E.: Parallel tempering algorithm for conformational studies of biological molecules. Chem. Phys. Let. 281, 140 - 150 (1997).

[10]: Hukushima, K. and Nemoto, K.: Exchange Monte Carlo method and application to spin glass simulations. J. Phys. Soc. (Jap.) 65, 1604-1608 (1996).

[11]: A. Kong, P. McCullagh, X.-L. Meng, D. Nicolae, and Z. Tan: A theory of statistical models for Monte Carlo integration. J. R. Statist. Soc. B 65, 585-618 (2003).

[12]: S. Kumar, D. Bouzida, R. H. Swendsen, P. A. Kollman and J. M. Rosenberg: The weighted histogram analysis method for free-energy calculations on biomolecules .1. the method. J. Comput. Chem. 13, 1011 - 1021 (1992).

[13]: Laio, A. and Parrinello, M.: Escaping free energy minima. Proc. Natl. Acad. Sci. USA 99, 12562 (2002).

[14]: E. Marinari and G. Parisi: Simulated Tempering: A New Monte Carlo Scheme. Euro. Phys. Let. 19, 451 (1992).

[15]: Antonia S. J. S. Mey, Hao Wu and Frank Noé: TRAM: Estimating equilibrium expectations from time-correlated simulation data at multiple thermodynamic states. http://arxiv.org/abs/1407.0138 (2014).

[16]: David D. L. Minh and John D. Chodera: Optimal estimators and asymptotic variances for nonequilibrium path-ensemble averages. J. Chem. Phys. 131, 134110 (2009).

[17]: J.-H. Prinz, J. D. Chodera, V. S. Pande, W. C. Swope, J. C. Smith and F. Noé: Optimal use of data in parallel tempering simulations for the construction of discrete-state Markov models of biomolecular dynamics. J. Chem. Phys. 134, 244108 (2011).

[18]: Edina Rosta and Gerhard Hummer: Free energies from dynamic weighted histogram analysis using unbiased Markov state model. J. Chem. Theo. Comput. 11, 276-285 (2015).

[19]: Michael R. Shirts and John D. Chodera: Statistically optimal analysis of samples from multiple equilibrium states. J. Chem. Phys. 129, 124105 (2008).

[20]: Marc Souaille and Beno{\^\i}t Roux: Extension to the weighted histogram analysis method: combining umbrella sampling with free energy calculations. Comput. Phys. Commun. 135, 40-57 (2001).

[21]: Sugita, Y. and Okamoto, Y.: Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 314, 141-151 (1999).

[22]: Z. Tan, E. Gallicchio, M. Lapelosa and R. M. Levy: Theory of binless multi-state free energy estimation with applications to protein-ligand binding. J. Chem. Phys. 136, 144102 (2012).

[23]: Torrie, G. M. and Valleau, J. P.: Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling. J. Comp. Phys. 23, 187-199 (1977).

[24]: H. Wu, A. S. J. S. Mey, E. Rosta and F. Noé: Statistically optimal analysis of state-discretized trajectory data from multiple thermodynamic states. J. Chem. Phys. 141, 214106 (2014).

[25]: Hao Wu and Frank Noé: Optimal estimation of free energies and stationary densities from multiple biased simulations. SIAM Multiscale Model. Simul. 12, 25-54 (2014).