% PRL look and style (easy on the eyes)
%\documentclass[aps,pre,nofootinbib,superscriptaddress,twocolumn,showkeys]{revtex4-1}
% Two-column style (for submission/review/editing)
\documentclass[aps,pre,preprint,nofootinbib,superscriptaddress,onecolumn]{revtex4-1}
% PRL look and style (easy on the eyes)
%\documentclass[aps,nofootinbib,superscriptaddress,twocolumn,showkeys]{revtex4-1}

\usepackage{palatino}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage[colorlinks=true,citecolor=blue,linkcolor=blue]{hyperref}

%\graphicspath{{fig/}}

\bibliographystyle{apsrevlong}

% italicized boldface for math (e.g. vectors)
\newcommand{\bfv}[1]{{\mbox{\boldmath{$#1$}}}}
% non-italicized boldface for math (e.g. matrices)
\newcommand{\bfm}[1]{{\bf #1}}          

\newcommand{\expect}[1]{\left \langle #1 \right \rangle}                % <.> for denoting expectations over realizations of an experiment or thermal averages

% vectors
\renewcommand{\vec}[1]{{\mbox{\boldmath{$#1$}}}}
\newcommand{\q}{\vec{q}} 				% coordinates
\newcommand{\p}{\vec{p}} 				% momenta
\newcommand{\z}{\vec{z}} 				% phase space point
\newcommand{\T}{\mathrm{T}}            % T used in matrix transpose

\newcommand{\qed}{\hfill \mbox{\raggedright \rule{.07in}{.1in}}}
\newcommand{\DeltaT}{\Delta t}

\DeclareMathOperator{\diag}{diag}
\def\argmax{\mathop{\text{argmax}}} 

% For red marked passages. Can later easily be switched of!
\newcommand{\markRed}[1] {{\color{red} #1}}

%% DOCUMENT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

%% TITLE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{Optimal use of data in parallel tempering simulations for the construction of discrete-state Markov models of biomolecular dynamics}

\author{Jan-Hendrik Prinz}
 \email{jan-hendrik.prinz@fu-berlin.de}
 \affiliation{Institute for Scientific Computing (IWR), University of Heidelberg, Im Neuenheimer Feld 368, 69126 Heidelberg, Germany and DFG Research Center Matheon, FU Berlin, Arnimallee 6, 14195 Berlin, Germany}

\author{John D. Chodera}
 \email{jchodera@berkeley.edu}
 \affiliation{Research Fellow, California Institute of Quantitative Biosciences (QB3), University of California, Berkeley, 260J Stanley Hall, Berkeley, California 94720, USA}

\author{Vijay S. Pande}
 \email{pande@stanford.edu}
 \affiliation{Department of Chemistry, Stanford University, Stanford, CA 94305}
 
\author{William C. Swope}
 \email{swope@almaden.ibm.com}
 \affiliation{IBM Almaden Research Center, 650 Harry Road, San Jose, California 95120}

\author{Jeremy C. Smith}
 \email{smithjc@ornl.gov}
 \affiliation{UT/ORNL Center for Molecular Biophysics, Oak Ridge National Laboratory P.O.Box 2008 Oak Ridge TN 37831-6164, USA}

\author{Frank No\'{e}}
 \email{noe@math.fu-berlin.de}
 \affiliation{DFG Research Center Matheon, FU Berlin, Arnimallee 6, 14195 Berlin, Germany}

\date{\today}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ABSTRACT/pacs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}
Parallel tempering (PT) molecular dynamics (MD) simulations have been extensively investigated as a means of efficient sampling of the configurations of biomolecular systems. Recent work has demonstrated how the short physical trajectories generated in PT simulations of biomolecules can be used to construct Markov models describing biomolecular dynamics at each simulated temperature. While this approach describes the temperature-dependent kinetics, it does not make optimal use of all available PT data, instead estimating the rates at a given temperature using only data from that temperature. This can be problematic, as some relevant transitions or states may not be sufficiently sampled at the temperature of interest, but might be readily sampled at nearby temperatures. Further, the comparison of temperature-dependent properties can suffer from the false assumption that data collected from different temperatures are uncorrelated. We propose here a strategy in which, by a simple modification of the PT protocol, the harvested trajectories can be reweighted, permitting data from all temperatures to contribute to the estimated kinetic model. The method reduces the statistical uncertainty in the kinetic model relative to the single temperature approach and provides estimates of transition probabilities even for transitions not observed at the temperature of interest. Further, the method allows the kinetics to be estimated at temperatures other than those at which simulations were run. We illustrate this method by applying it to generating a Markov model of the conformational dynamics of the solvated terminally-blocked alanine peptide.

\emph{Keywords: temperature reweighting; temperature-dependent dynamics; parallel tempering; molecular dynamics; markov state model.}
\end{abstract}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INTRODUCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{section:introduction}

Biological macromolecules are not static structures, but are driven by thermal motion and interactions with their molecular environment, undergoing conformational fluctuations and changing conformational states. The characterization of the statistical conformational dynamics of biomolecules is essential to understanding how these molecules work as molecular machines.

Often, a separation of timescales of characteristic dynamical relaxation times gives rise to the existence of metastable conformational states, such that the biomolecule remains in any one of these states for a long time before making a rapid transition to another state.
A wealth of experimental data now supports the existence of such states, including NMR~\cite{kern:science:2002,kern:nature:2005,Feng:2005p6514}, fluorescence emission~\cite{Min:2005p4264,Smith:2008p4287}, energy transfer~\cite{zhuang:2003a,Matagne:1997p6517}, correlation spectroscopy~\cite{Neuweiler:2005p5935,Goldbeck:1999p6516}, and non-equilibrium perturbation experiments~\cite{Smith:2008p4287}. 
Developing a quantitative understanding of what gives rise to these conformational states and the interactions that govern transitions between them will have a significant impact on our understanding of many biological processes, such as, for example, signaling events, enzyme regulation, allostery, and drug design with conformationally flexible molecules.

Sampling the underlying phase space by straightforward molecular dynamics simulation often suffers from the problem that the timescales of conformational changes can be orders of magnitude larger than simulation times accessible using current computational resources. PT molecular dynamics simulation has been an effective and thus popular approach to overcoming the issue of convergence in molecular simluations, by allowing replicas to heat up and overcome enthalpic barriers as the simulation proceeds while still sampling from an appropriate equilibrium distribution~\cite{geyer-thompson:j-am-stat-assoc:1995:parallel-tempering,hukushima-nemoto:j-phys-soc-jpn:1996:parallel-tempering,hansmann:chem-phys-lett:1997:parallel-tempering-monte-carlo,Sugita_CPL314_141}. At the same time this approach permits an analysis of the temperature dependence of properties of interest, which is especially important for comparisons with certain experimental results (e.g.~melting curves, heat capacities)~\cite{Jager:2001p6248}. Although PT molecular dynamics produces unphysical replica trajectories, the short physical trajectories in between the exchanges can provide useful dynamical information. If the PT simulation is well-equilibrated, these initial configurations of the short trajectory segments will be sampled from the equilibrium at their corresponding temperatures.

Markov models provide a way of modeling the slow conformational dynamics of biomolecules based on short simulations~\cite{Buchete:2008p2351,Noe:2007p867,Chodera:2006p2352,Muff:2009p6071,Chodera:2007p2463,Hinrichs:2007p5942,Singhal:2005p2464,Swope:2004p907,Swope:2004p1280}.
In these models conformational states are envisioned as disjoint but connected regions of configuration space. The biomolecule spends long times within individual regions before undergoing rapid stochastic transitions between them. 
If a separation of timescales exists between fast relaxation times within and slow equilibration between regions, the inter-state dynamics can be well described by a Markov model in discrete timesteps $\tau$, where coarse graining in time is required as the discretization in space prohibits the characterization of relaxation processes faster than $\tau$. If the system is partitioned into its metastable states, $\tau$ is related to the time required to overcome internal barriers within each conformational state. However, it has been recently shown that, even in the absence of many metastable states, a Markovian model can well-approximate the dynamics at long times, with this approximation error decreasing as the number of states is increased~\cite{sarich:errorindiscretization}.

Recently, Buchete and Hummer have shown that both thermodynamic and kinetic properties can be estimated over the range of temperatures by constructing Markov models using the short physical trajectories generated from PT simulations~\cite{Buchete:2008p2354}.
However, if a complete description of dynamics across the entire thermodynamically relevant configurational space at a given temperature is desired, one quickly runs into problems if use is made of trajectories only from the temperature of interest, as some states that are sampled at other temperatures may not be well-sampled at the single temperature~\cite{Buchete:2008p2354}. One would like to make use of the data collected at \emph{all} temperatures to characterize the kinetic behavior in all regions sampled over the full range of temperatures spanned by the PT simulation in a manner similar to equilibrium reweighting~\cite{Ferrenberg:1989p3671,bartels:1997a,Kumar:1992p860,gallicchio:2005a,Bartels:2000p4835,Shirts:2008p2350}.

Here, we propose a method for integrating MD data from all temperatures by making use of \emph{dynamical reweighting}~\cite{Chodera:2009p6073}, allowing a smooth, continuous, and differentiable estimate of the transition probabilities at any temperature without requiring the assumption of any kinetic model (such as Arrhenius kinetics~\cite{Muff:2009p6071}) and taking advantage of the increased transition rates at higher (or, for transitions with entropic barriers, lower) temperatures. Reweighting methods (such as histogram-based~\cite{Ferrenberg:1989p3671,bartels:1997a,Kumar:1992p860,gallicchio:2005a} or histogram-free~\cite{Bartels:2000p4835,Shirts:2008p2350} approaches) allow the use of samples collected from multiple distributions to provide an improved estimate of the expectation value of some static property at the distribution of interest, and have been used extensively in the analysis of equilibrium thermodynamic properties in replica-exchange simulations~\cite{Mitsutake:2001p3675}.

Dynamical reweighting has recently been proposed as a way of estimating dynamical properties (such as correlation functions) using an asymptotically optimal estimator, simultaneously providing a good estimate of the statistical error~\cite{Chodera:2009p6073}. Here, we show how dynamical reweighting can be used to estimate transition probabilities (and their statistical uncertainties) for the construction of a Markov model as a smooth function of temperature, making use of data from all temperatures. This has the advantage of producing a useful Markov model at \emph{any} temperature containing the dependence of kinetic properties on temperature, and providing an assessment of the error in the model.

We illustrate this approach for the standard test case of the terminally-blocked alanine peptide in explicit solvent. A Markov model constructed from short (6 ps) trajectories from each state has been previously shown to accurately describe the kinetics of this system at 302~K~\cite{Chodera:2006p2352}. This peptide system presents a challenge for estimators based on individual temperatures, due to the presence of highly metastable states with very high free energies relative to the most populated states. These states are poorly sampled at temperatures near 300 K, even though their temporal behavior can dominate the non-equilibrium relaxation kinetics at this temperature. Finally, we determine whether using all the data using reweighting produces substantially improved kinetic models at this particular temperature and across the full range of temperatures in the PT simulation.

This paper is organized as follows:
In Section~\ref{section:theory} we review the theory behind Markov models of multistate conformational dynamics. We then show in Section~\ref{section:dynamical-reweighting} how dynamical reweighting can be used to estimate temperature-dependent transition probabilities and rates for a given state decomposition. Finally, we illustrate the method in Section~\ref{section:application-to-alanine-dipeptide} by applying it to a six-state decomposition of the terminally-blocked alanine peptide, and compare the results to the approach of Buchete and Hummer~\cite{Buchete:2008p2354}, in which Bayesian estimates of transition probabilities are obtained from individual temperatures extracted from the PT simulation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THEORY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\color{red}
\section{Theory}
\label{section:theory}
\subsection{Markov Models}
Consider a system that evolves according to some stationary dynamical process. Let $\Omega$ be the configuration space with a complete decomposition~\cite{Swope:2004p907} into $M$ disjoint sets $\Gamma_{i} \subset \Omega$, $i \in \{ 1,\ldots,M \}$, such that
\begin{equation}
	\bigcup\limits_{i=1}^{M} \Gamma_{i} = \Omega; \; \Gamma_{i} \cap \Gamma_{j} = \emptyset, \; i \neq j.
\end{equation}
For convenience we also define indicator functions $\chi_{i}(\bfv{q})$ for points in configuration space $\bfv{q}\in\Omega$, such that
\begin{equation}
	\chi_{i}(\bfv{q}) \equiv 
	\begin{cases}
		1 & \text{if} \; \bfv{q} \in \Gamma_{i} \\
		0 & \text{otherwise}
	\end{cases} ,
\end{equation}
i.e., the function assumes the value of unity if $\bfv{q}$ belongs to set $\Gamma_{i}$, and zero otherwise. 

Based on this discretization of configuration space, we can define a $M\times M$ row-stochastic transition matrix $\bfm{T}(\tau)$ with conditional probabilities of finding the system in state $j$ a time $\tau$ after it was originally in state $i$,
\begin{eqnarray}
T_{ij}(\tau) & \equiv & \mathbb{P}\left(\bfv{q}(\tau) \in \Gamma_{j} \mid \bfv{q}(0) \in \Gamma_{i}\right) \\
 & = & \frac{\left\langle \chi_{i}(0) \chi_{j}(\tau) \right\rangle}{\left\langle \chi_{i} \right\rangle}, \label{equation:transition-matrix-definition}
\end{eqnarray}
and introduce the notation $\chi_{i}(t) \equiv \chi_{i}(\bfv{q}(t))$, where the dynamics is assumed to be governed by a stationary (time-independent) process.
We aim here to construct a discrete-time, discrete-space Markov model that approximates the long-time dynamics of the system by virtue of
\begin{equation}
	\bfv{p}(t + k \tau) \approx \bfv{p}(t) \bfm{T}^k(\tau)
	\label{eq:approximateProbabilityTransport}
\end{equation}
with $\bfv{p}$ being the projection of some continuous distribution $\rho(\bfv{q})$ onto the discrete subsets $\Gamma_{i}$.

Eq.~\ref{eq:approximateProbabilityTransport} is only an approximation to the real dynamics due to the introduction of the coarse-graining of configuration space into $M$ discrete sets~\cite{Swope:2004p907,Swope:2004p1280}. 
Despite this, recent theoretical work has demonstrated that the approximation error introduced by the discretization can be made arbitrarily small by increasing the number of states $M$, increasing the lagtime $\tau$, or both~\cite{sarich:errorindiscretization}.
As a result, we can ensure the true dynamics is reproduced by a discrete-time, discrete-space Markov model to arbitrary desired precision through careful choice of $M$ and the sets $\Gamma_i$.
Although this is a crucial step in the construction of a Markov model of dynamics, the process of finding an optimal decomposition of state space and appropriate lagtime $\tau$ is beyond the scope of this paper, and has been extensively discussed elsewhere~\cite{Noe:2007p867,Chodera:2007p2463}.
\color{black}

\subsection{Estimating transition probabilities \emph{vs} estimating transition rates}
\label{section:probabilities-vs-rates}
\color{red}
Coarse-grained dynamics is equivalently described by the discrete-state, continuous-time \emph{master equation}
\begin{equation}
	\dot{\bfv{p}}(t) \approx \bfv{p}(t) \; \bfm{K} 
	\label{equation:master-equation}
\end{equation}
where $\bfm{K}$ denotes the $M \times M$ phenomenological rate matrix. For $i\neq j$, $K_{ij} > 0$ is the rate constant associated with the transition $i\rightarrow j$, while the diagonal elements are set to $K_{ii} = - \sum_{j\neq i}K_{ij}$ for conservation of probability~\cite{Buchete:2008p2354,Voter:1985p2344,Adams:1981p2340}. It is noteworthy that, although Eq.~\ref{equation:master-equation} is defined for all times $t \ge 0$, the time-evolution of $\bfv{p}(t)$ is accurately modeled only for times $t$ larger than the same Markov time $\tau$ appearing in Eq.~\ref{eq:approximateProbabilityTransport}.

While $\bfm{T}(\tau)$ can be straightforwardly estimated from a trajectory using Eq.~\ref{equation:transition-matrix-definition}, $\bfm{K}$ cannot, because inversion of the matrix exponential in equation
\begin{equation}
	\bfm{T}(\tau) = \exp( \bfm{K} \tau )
\end{equation}
is not unique for stochastic matrices unless $\bfm{T}(\tau)$ is also positive definite and reversible, and even in these cases, numerical issues can frustrate the computation of the matrix logarithm of a positive-definite reversible $\bfm{T}(\tau)$. One potential solution to these issues is to use Bayesian inference to estimate the likely rate matrix given data~\cite{Buchete:2008p2351}---we discuss this issue further in Section~\ref{section:reversible-rate-matrices}. For simplicity, we focus here on the estimation of temporally discrete transition probabilities $\bfm{T}(\tau)$ in the next few sections. 
\color{black}

\subsection{Estimation from temporally discrete trajectories}
\label{section:discrete-trajectories}
\markRed{
Consider a configuration-space trajectory $\bfv{q}(t) \in \Omega$, $t \in [0,t_\mathrm{max}]$, sampled at time intervals of $\DeltaT$.
We construct a temporally discrete version of the trajectory, $\bfv{q}_n \equiv \bfv{q}(n \Delta t)$, and introduce the discrete Markovian lagtime, $s \equiv \tau / \DeltaT \in \mathbb{N}$, and discrete trajectory length, $L \equiv t_\mathrm{max} / \DeltaT \in \mathbb{N}$.

If the lagtime $\tau$ is long enough, the statistical dynamics over times $\tau$ and longer can be well-approximated by a Markov chain \cite{Chodera:2006p6224,Swope:2004p1280,sarich:errorindiscretization}, and can we define the count matrix $\bfm{B}(\tau)$ by 
\begin{eqnarray}
	B_{ij}(\tau) & \equiv & \frac{1}{s}\sum\limits_{n = 0}^{L-s} \chi_{i}\left(\bfv{q}_{n}\right) \chi_{j}\left(\bfv{q}_{n + s}\right)
	\label{equation:count-matrix}
\end{eqnarray}
as the effective number of (potentially fractional) transitions observed among the states for a fixed lagtime $\tau$.
The associated likelihood $\mathcal{L}$ that a given transition matrix $\bfm{T}$ produces the observations stored in the count matrix $\bfm{B}$ is given by the multinomial distribution
\begin{equation}
	\mathcal{L}(\bfm{T}) = \mathbb{P}(\bfm{B} | \bfm{T}) \propto \prod \limits_{i,j = 1}^{M} T_{ij}^{B_{ij}}.
	\label{equation:likelihood}
\end{equation}
where, as a representative, we choose the unique transition matrix $\hat{\bfm{T}}(\tau)$ which maximizes this likelihood:
\begin{equation}
	\hat{T}_{ij}(\tau) = \left[\argmax_{\bfm{T}} \; \mathcal{L}(\bfm{T})\right]_{ij} = \frac{B_{ij}(\tau)}{\sum\limits_{k} B_{ik}(\tau)}.
\end{equation}
Alternatively, we can use the state-to-state time-correlation function $C_{ij}(\tau)$~\cite{Swope:2004p1280,Swope:2004p907} given by
\begin{eqnarray}
	C_{ij}(\tau) &\equiv& \expect{\chi_i(0) \, \chi_j(\tau)}
\end{eqnarray}
which can be estimated by $\hat{C}_{ij}(\tau)$ in a similar fashion 
\begin{eqnarray}
	\hat{C}_{ij}(\tau) & = & \frac{1}{L-s}\sum\limits_{n = 0}^{L-s} \chi_{i}(\bfv{q}_{n}) \chi_{j}(\bfv{q}_{n + s}) = \frac{s}{L-s} B_{ij}(\tau).
	\label{equation:definition-of-correlation}
\end{eqnarray}
Although dynamical reweighting can be formulated for different dynamical models~\cite{Chodera:2009p6073}, the present approach is based on Hamiltonian dynamics in the canonical ensemble, which is time-reversible, and thus equilibrium molecular dynamics fulfills detailed balance in state space. Consequently, for trajectories sampled from equilibrium, the correlation matrix will have a symmetric form $C_{ij}(\tau) = C_{ji}(\tau)$. In this case, we can use the estimator
\begin{equation}
	\hat{C}_{ij} = \frac{s}{2 (L-s)} (B_{ij} + B_{ji}) = \hat{C}_{ji}
	\label{eq:correlation-from-counts}
\end{equation}
and write the transition matrix estimate $\hat{\bfm{T}}(\tau)$ in terms of the correlation matrix estimate $\hat{\bfm{C}}(\tau)$ as
\begin{eqnarray}
 \hat{T}_{ij}(\tau) & = & \frac{\hat{C}_{ij}(\tau)}{\sum\limits_{k=1}^M \hat{C}_{ik}(\tau)}
 \label{eq:transitionmatrix-from-correlation}
\end{eqnarray}
which will also fulfill detailed balance~\cite{Swope:2004p1280,Swope:2004p907}.
}
\subsection{Transition probabilities from dynamical reweighting}
\label{section:dynamical-reweighting}

\markRed{
We now demonstrate how transitions observed at all temperatures can be used to infer transition probabilities at any temperature of interest through the use of a reweighting procedure. 
We consider the specific case of a system in the canonical ensemble at inverse temperature $\beta \equiv (k_{B} T)^{-1}$, where dynamics from a canonical distribution of initial phase space points obeys Hamiltonian dynamics over the observed lagtime $\tau$.
This corresponds to the case of ``canonical distribution of Hamiltonian trajectories'' examined in detail in Ref.~\cite{Chodera:2009p6073}.

With this definition of dynamics, the (now temperature-dependent) state-to-state correlation functions $C_{ij}(\tau;\beta)$ can be expressed as Boltzmann-weighted expectation functions.
Defining $\bfv{z} \equiv (\bfv{q},\bfv{p})$ as a point in phase space, we can write,
\begin{eqnarray}
	C_{ij}(\tau;\beta) = \frac{1}{Z(\beta)} \int d \bfv{z} \, e^{ -\beta \mathcal{H}(\bfv{z}_{0}) } \chi_i(\bfv{q}_{0}) \chi_j(\bfv{q}_{\tau}) 
	\label{equation:temperature-dependent-correlation-function},
\end{eqnarray}
where $Z(\beta)$ is the complete partition function of both kinetic and potential energies and $\bfv{z}_\tau \equiv (\bfv{q}_\tau, \bfv{p}_\tau)$ is the phase space point produced by evolving Hamiltonian dynamics a time $\tau$ from an initial phase space point $\bfv{z}_0 \equiv (\bfv{q}_0, \bfv{p}_0)$.

Suppose we have a set of $N$ Hamiltonian phase space trajectories $\z_{n}(t)$, $n \in \{ 1,\ldots,N \}$ of length $t_{\text{max}}$, $t \in [0,t_{\text{max}}]$ sampled from equilibrium at $K$ temperatures $\beta_k$, $k  \in \{ 1,\ldots,K \}$.
For convenience, we group the $N$ trajectories into subsets $\mathcal{Q}_{k} \subset \{1,\ldots,N\}$ according to which temperature $\beta_k$ their initial phase space points $\z_{n}(0)$ were drawn from.
For the reweighting procedure, only the number of trajectories $N_{k} \equiv \left| \mathcal{Q}_{k} \right| $ sampled from the temperature $\beta_{k}$ is relevant, and not the direct association of a specific trajectory $\bfv{z}_{n}$ with the temperature $\beta_k$ it was sampled from~\cite{Shirts:2008p2350,Chodera:2009p6073}.

By the application of \emph{dynamical reweighting}~\cite{Chodera:2009p6073}, a correlation function $C_{ij}(\tau;\beta)$ can then be estimated using the entire set of trajectories at all temperatures with
\begin{eqnarray}
	\hat{C}_{ij}(\tau; \beta) & \approx & \sum\limits_{n=1}^{N} w_{n}(\beta) \cdot \hat{C}_{ij}^{(n)}(\tau;\beta)
	\label{correlation-matrix-from-DR}
\end{eqnarray}
where the individual trajectory contributions $\hat{C}_{ij}^{(n)}(\tau;\beta)$ are estimated using Eq.~\ref{eq:correlation-from-counts} and the normalized trajectory weights $w_{n}(\beta)$ can be computed by
\begin{eqnarray}
	w_{n}(\beta) &=& \hat{Z}(\beta)^{-1} \left[ \sum_{k=1}^K N_{k} \hat{Z}_{k}^{-1} e^{ - (\beta_{k} - \beta) E_{n} } \right]^{-1} 
	\label{equation:trajectory-weights},
\end{eqnarray}
where the normalization constants $\hat{Z}(\beta)$ are 
chosen that
\begin{equation}
	\sum_{n=1}^{N} w_{n}(\beta) = 1
\end{equation}
holds
%determined by
%\begin{eqnarray}
%	\hat{Z}(\beta) &\equiv& \sum_{n=1}^N \left[ \sum_{k=1}^K N_{k} \hat{Z}_{k}^{-1} e^{ - (\beta_{k} - \beta) E_{n} } \right]^{-1}
%\end{eqnarray}
and $E_{n} \equiv H(\z_{n}(0))$ denotes the total energy of the system in trajectory $\z_{n}$, which is constant over trajectories for Hamiltonian dynamics. 

The constants $\hat{Z}_{k}$ need to be determined from the solution of a set of $K$~self-consistent equations indexed by $ i \in \{ 1, \ldots, K \}$
\begin{eqnarray}
	\hat{Z}_{i} & = & \sum_{n=1}^N \left[ \sum_{k=1}^K N_{k} \hat{Z}_{k}^{-1} e^{ - (\beta_{k} - \beta_i) E_{n} } \right]^{-1}
	\label{eq:free-energies}
\end{eqnarray}
which can be obtained efficiently in a number of ways~(see Appendix~\ref{appendix:solving-self-consistent-equations}), although it is often necessary to work with logarithmic representations to avoid numerical instability. The choice of weights $w_{n}$ in Eq.~\ref{equation:trajectory-weights} gives an asymptotically optimal (i.e.~lowest variance) estimate of the
correlation function in Eq.~\ref{equation:temperature-dependent-correlation-function}. 
A detailed exposition of dynamical reweighting for the estimation of correlation functions is presented in Ref.~\cite{Chodera:2009p6073}.

Finally, the row-stochastic transition matrix estimate $\hat{\bfm{T}}(\tau;\beta)$ is computed from Eq.~\ref{correlation-matrix-from-DR} using Eq.~\ref{eq:transitionmatrix-from-correlation}, where the symmetry of $\hat{\bfm{C}}(\tau;\beta)$ results in a reversible transition matrix estimate $\hat{\bfm{T}}(\tau;\beta)$ (i.e., it will satisfy detailed balance).
}

\subsection{Estimation of uncertainties in transition probabilities}
\label{appendix:propagation-of-error}
\markRed{
For a given temperature $\beta$ and lagtime $\tau$, the statistical uncertainty in $\hat{C}_{ab} \equiv \hat{C}_{ab}(\tau;\beta)$ can be estimated in a straightforward manner~\cite{Shirts:2008p2350,Chodera:2009p6073}. We start with the $N\times K$ weight matrix $\bfm{W}$, the elements of which are given by
\begin{equation}
	w_{nk}  \equiv  w_{n}(\beta_{k})
\end{equation}
where $k \in \{ 1,\ldots,K \}$ runs over the set of temperatures $\beta_k$ and $n \in \{ 1,\ldots,N \}$ over all trajectories.
Augmenting $\bfm{W}$ by three additional columns, indexed by the letters $\text{A}$, $\text{X}$, and $\text{Y}$, defined as
\begin{eqnarray}
	w_{n\text{A}} & = & w_n(\beta) \\
	w_{n\text{X}} = w_{nA} \frac{\hat{C}_{ab}^{(n)}}  {\hat{C}_{ab}} & ; & 	w_{n\text{Y}} = w_{nA} \frac{\hat{C}_{a'b'}^{(n)}}{\hat{C}_{a'b'}} ,
\end{eqnarray}
the covariance of two correlation matrix elements $\hat{C}_{ab}$ and $\hat{C}_{a'b'}$ can be estimated by~\cite{Chodera:2009p6073}
\begin{equation}
\delta \hat{C}_{ab} \delta \hat{C}_{a'b'} \approx \hat{C}_{ab} \, \hat{C}_{a'b'} \left[ \hat{\Theta}_{AA} - \hat{\Theta}_{AY} - \hat{\Theta}_{XA} +\hat{\Theta}_{XY} \right] \label{equation:covariance-of-estimates}
\end{equation}
where the covariance matrix estimate $\hat{\bfm{\Theta}}$ is computed as~\cite{Chodera:2009p6073}
\begin{eqnarray}
\hat{\bfm{\Theta}} &\equiv& \bfm{W}^T [\bfm{I}_N - \bfm{W} \bfm{N} \bfm{W}^\T]^+ \bfm{W}, \label{equation:asymptotic-covariance}
\end{eqnarray}
with $\bfm{I}_{N}$ being the identity matrix of rank $N$, $\bfm{N}=\mathrm{diag}(N_1, \ldots, N_K, 0, 0, 0)$ and the generalized inverse denoted by $[]{}^{+}$.

The uncertainty in the estimated transition probabilities $\delta^{2}\hat{T}_{ij}$ can be approximated as a function of the uncertainty in the correlation matrix elements, $\delta \hat{C}_{ab} \delta \hat{C}_{a'b'}$, by a first-order Taylor expansion about the mean $\expect{\hat{T}_{ij}}$,
\begin{eqnarray}
	\delta^2 \hat{T}_{ij} & \equiv & \expect{\left(\hat{T}_{ij} - \expect{\hat{T}_{ij}}\right)^2} \nonumber \\
&\approx& \sum_{a,a^{'},b,b^{'}=1}^M \left[ \frac{\partial \hat{T}_{ij}}{\partial \hat{C}_{ab}} \right] \, \left[ \frac{\partial \hat{T}_{ij}}{\partial \hat{C}_{a'b'}} \right] \, \delta \hat{C}_{ab} \delta \hat{C}_{a'b'}.
\label{equation:uncertainty-estimation-main}
\end{eqnarray}
Using Eq.~\ref{equation:transition-matrix-definition}, the sensitivity of $\hat{T}_{ij}$ to the correlation matrix element $\hat{C}_{ab}$ is given by
\begin{equation}
\frac{\partial \hat{T}_{ij}}{\partial \hat{C}_{ab}} = \frac{\delta _{aj} \delta _{bi}+\delta _{ai} \delta _{bj}-\delta
   _{ab}}{\hat{C}_i} - \frac{\hat{C}_{ij} \left(-M \delta _{ab}+\delta _{ai}+\delta
   _{bi}\right)}{\hat{C}_i^2},
   \label{equation:estimation-of-transition-matrix-elements}
\end{equation}
with $\hat{C_{i}} \equiv \sum_{j=1}^M \hat{C}_{ij}$ representing the row sums of the correlation matrix $\hat{\bfm{C}}$.

The complete expression for the variance in the transition probabilities $\delta^2 \hat{T}_{ij}$ in Eq.~\ref{equation:uncertainty-estimation-main} can be evaluated using Eqs.~\ref{equation:covariance-of-estimates},~\ref{equation:asymptotic-covariance} and~\ref{equation:estimation-of-transition-matrix-elements}. The final result is complex but still calculable, as demonstrated in section~\ref{section:application-to-alanine-dipeptide}. 
A detailed description of the procedure for computing statistical uncertainties for arbitrary properties can be found in Ref.~\cite{Chodera:2009p6073,Shirts:2008p2350}.
}
\subsection{Modified parallel tempering protocol}
\label{section:modified-parallel-tempering-protocol}
\markRed{
We employ a modified parallel tempering (PT) protocol in which a set of $K\times L$ Hamiltonian trajectory segments $\z_{kl}(t)$ of uniform length $t_{\text{max}}$ is generated~\cite{Chodera:2009p6073} with temperature index $k \in \{ 1,\ldots,K \}$ and iteration $l \in \{ 1, \ldots, L \}$. For the Markov property to hold in the later estimation, the length $t_{\text{max}}$ between proposed exchanges in the PT protocol needs to exceed the lagtime $t_{\text{max}} \ge \tau$. 
We start by assuming that some process was used to generate the initial phase space points $\z_{k0}(0)$ for the first set of  trajectories ($l=0$) from equilibrium within the canonical (NVT) ensemble at each corresponding inverse temperature $\beta_k$
\begin{eqnarray}
	\mathbb{P}\left(\z_{k0}(0) ; \beta_k\right) & = & \frac{1}{Z(\beta_k)} \, e^{-\beta_{k} \mathcal{H}(\z_{k0}(0))}.
\end{eqnarray}
These initial phase space points may be obtained, for example, by a standard PT protocol, or by running the modified protocol for a number of iterations starting from one or more arbitrary initial configurations. 

In each subsequent iteration ($l > 0$), Hamilton's equations of motion are used to propagate all replicas using a symplectic integrator with sufficiently small timesteps to generate trajectories $\z_{kl}(t)$ of length $t_{\text{max}}$. At $t = t_{\text{max}}$, we propose exchanges between the final configurations $\z_{il}(t_{\text{max}})$ and $\z_{jl}(t_{\text{max}})$ of neighboring temperatures $\beta_i$ and $\beta_{i \pm 1}$, starting from the highest temperature down to the lowest one in odd iterations and in reverse order in even ones~\cite{Chodera:2006p2352}.\footnote{Note that other exchange proposal schemes can be used, provided the resulting algorithm satisfies the condition of ``balance'' (which is less strict than detailed balance)~\cite{Earl:2005p2409}.} 
The Metropolis-like probability~\cite{Sugita_CPL314_141} for accepting or rejecting the exchange $P_\mathrm{exch}$ depends on the potential energies of the final configurations $U_{k} \equiv U(\z_{kl}(t_{\text{max}}))$ with
\begin{eqnarray}
	P_\mathrm{exch}(U_i, \beta_i; U_j, \beta_j) &=& \min \left\{ 1, e^{-(\beta_i - \beta_j) (U_j - U_i)} \right\}.
\end{eqnarray}
}
\color{red}
Regardless of whether the exchange is accepted or rejected, we reassign the velocities according to the Maxwell-Boltzmann distribution~\cite{andersen:1980:jcp:andersen-thermostat} at the new (or old, if rejected) temperatures, and denote the new phase space points from which the next iteration can be carried out as $\z_{k(n+1)}(0)$ (see proof in Appendix~\ref{appendix:parallel-tempering-stationary-distribution}).
This satisfies the conditions defined by Okamoto~\cite{Sugita_CPL314_141} in order for the kinetic energies to not appear in $P_\mathrm{exch}$ and is equivalent to rescaling the velocities for accepted exchanges and then applying a massive collision for the Andersen thermostat~\cite{andersen:1980:jcp:andersen-thermostat}. 
The reason for reassignment of velocities instead of rescaling is that when using Hamiltonian trajectories to propagate dynamics, no thermostatting would otherwise take place. 
Without this velocity reassignment, the use of Hamiltonian trajectories (even if velocity rescaling were performed after exchanges) would mean that no or minimal thermostatting would take place, generating an improper ensemble; velocity reassignment ensures the canonical ensemble is generated.
\color{black}

\subsection{Bayesian estimation of transition probabilities from a single temperature}
We also consider two Bayesian methods for estimation of the transition matrices and rate matrices using data collected from a single temperature. Both methods sample transition probabilities or rates according to the same likelihood function in Eq.~\ref{equation:likelihood}, but employ different model parameterizations and, more importantly, different prior probability distributions.

\subsubsection{Reversible Transition Matrices}
\label{section:transition-matrix-estimation}
We first consider the approach described in Ref.~\cite{Noe:2008p2459} to infer transition matrices that satisfy detailed balance. 
Starting with an observation represented by the fractional count matrix $\bfm{B}$, the posterior probability a transition matrix $\bfm{T}$ was responsible for generating this observation is given by
\begin{equation}
\mathbb{P}(\bfm{T}|\bfm{B}) \propto \mathbb{P}(\bfm{B}|\bfm{T}) \; \mathbb{P}(\bfm{T}) =  \prod_{i,j = 1}^{M} T_{ij}^{B_{ij}} \; \mathbb{P}(\bfm{T}).
\label{eq:frank-likelihood}
\end{equation}
As the prior, $\mathbb{P}(\bfm{T})$, we choose a Dirichlet distribution for each row which adds no additional observations to the likelihood probability
\begin{equation}
	\mathbb{P}(\bfm{T}) \equiv \prod_{i,j=1}^{M} T_{ij}^{-1}.
\end{equation}
Furthermore, we restrict ourselves to transition matrices that fulfill detailed balance, i.e., are reversible with respect to the stationary distribution $\bfm{\pi}$:
\begin{equation}
	\pi_{i} T_{ij} = \pi_{j} T_{ji}
\end{equation}
Here, the distribution in Eq.~\ref{eq:frank-likelihood} was sampled using a Markov chain Monte Carlo procedure described in Ref.~\cite{Noe:2008p2459}.

\subsubsection{Reversible Rate Matrices}
\label{section:reversible-rate-matrices}

\color{red}
We also consider a second approach to estimating transition probabilities from individual temperatures using Bayesian estimation.
Here, we sample rate matrices $\bfm{K}$ (discussed in section~\ref{section:probabilities-vs-rates}) with elements $K_{ij}>0$ for $i\neq j$ and $K_{ii} = - \sum_{j\neq i} K_{ij}$ using the approach described in Ref.~\cite{Buchete:2008p2354}.
This approach does not estimate the transition probabilities directly, but instead uses a parametric form of a reversible rate matrix $\bfm{K}$, that uses the logarithms of the elements in the upper-right triangular matrix $K_{ij}$ for $j > i$ (without diagonal entries) and the equilibrium distribution $\pi_{i}$, thus assuring a rate matrix with non-positive eigenvalues and also positive off-diagonal rates.
\color{black}

The posterior in Eq.~\ref{eq:frank-likelihood}, written in terms of the rate matrix $\bfm{K}$, is given by
\begin{equation}
	\mathbb{P}(\bfm{K}|\bfm{B}) \propto \mathbb{P}(\bfm{B}|\bfm{K}) \mathbb{P}(\bfm{K}) = \prod_{i,j = 1}^{M} \exp \left( \tau \bfm{K}\right) _{ij}^{B_{ij}} \mathbb{P}(\bfm{K})
	\label{eq:hummer-likelihood}
\end{equation}
where the prior is uniform in $\ln K_{ij}, j>i$ and $\ln \pi_{i}$ and the detailed balance constraint
\begin{equation}
	\pi_{i} K_{ij} = \pi_{j} K_{ji}
\end{equation}
holds.
After sampling rate matrices with a Metropolis Monte Carlo scheme~\cite{Buchete:2008p2354} the related set of transition matrices with the lagtime $\tau$ is computed by
\begin{equation}
	\bfm{T}(\tau) = \exp \left( \tau \bfm{K}\right)
	\label{eq:ratematrix-to-transitionmatrix}
\end{equation}
All methods with their abbreviations and colors used consistently in the text and figures in this article are listed in Table~\ref{table:methods-colors}.
\begin{table}
\begin{ruledtabular}
\begin{tabular}{lccc}
Method & Abbreviation & Color \\ \hline
Transition matrix estimation~\cite{Noe:2008p2459} & [TE] & RED\\
Rate matrix estimation~\cite{Buchete:2008p2354} & [RE] & GREEN\\
Dynamical reweighting~\cite{Chodera:2009p6073}& [DR] & BLUE\\
Shooting trajectories at 302 K~\cite{Chodera:2006p2352} & [ST] & BLACK 
\end{tabular}
\end{ruledtabular}
\caption{\label{table:methods-colors}Table of methods used for transition probability or rate estimation with their corresponding abbreviations and colors used consistently throughout this paper.}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% APPLICATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Application to terminally-blocked alanine dipeptide}
\label{section:application-to-alanine-dipeptide}
\subsection{System Setup}

\markRed{
To illustrate the construction of temperature-dependent Markov models using dynamical reweighting, we estimated the transition probabilities between conformational states for the terminally-blocked alanine peptide (Ace-Ala-Nme) (see Fig.~\ref{fig:energyLandscape}) in explicit solvent from a parallel tempering (PT) molecular dynamics simulation. 
This dataset\footnote{The alanine dipeptide parallel tempering and kinetics datasets are available online for download at \url{https://simtk.org/home/alanine-dipeptide/}.} was published previously as part of a study that explored the suitability of a Markov model for describing dynamics at 302 K~\cite{Chodera:2006p2352}; here, we make use of this dataset to facilitate comparison to an independent trajectory set that appears in that reference.
The PT dataset consists of an ensemble of 501 Hamiltonian trajectories 20 ps in length at each of 40 temperatures, spanning 273 to 600 K, with peptide configurations stored every 0.1 ps. 
The temperatures were exponentially spaced to ensure good overlap in the potential and total energy distributions between neighboring temperatures and reasonably high exchange probabilities (see Fig.~\ref{figure:energydistribution}).
}

A leapfrog Verlet integrator~\cite{verlet:1967:phys-rev:verlet-integrator-1,verlet:1967:phys-rev:verlet-integrator-2,berendsen:1977:j-comput-phys:SHAKE} (with bonds involving hydrogen atoms constrained) was used to produce the dynamical trajectories. The fluctuation in total energy averaged over all 20 ps trajectories at each temperature was minimal and the drift negligible (see Table~\ref{table:pt-energy-rms-and-drift}). The production run followed a 1 ns equilibration phase during which exchanges were attempted at 1 ps intervals, ensuring that all initial configurations were drawn from equilibrium at their respective temperatures.
\begin{figure*}[]
\includegraphics[width=\textwidth]{fig/plot_EnergyOverlap_InitialEnergy_Gauss.pdf}
\caption{\label{figure:energydistribution} {\bf Distribution of trajectory total energies from parallel tempering simulation.} 
Distributions for trajectories sampled at temperatures 302~K, 425~K, and 600~K---for which single-temperature Bayesian analysis convergence properties are shown in Fig.~\ref{fig:convergence-of-bayesian-methods}---are highlighted.}
\end{figure*}
\begin{table}
\begin{ruledtabular}
\begin{tabular}{ccc}
Temp [K] & SD [kcal/mol] & Drift [kcal/(ps mol)] \\ \hline
302 & $0.214 \pm 0.014$  & $0.0056 \pm 0.0006 $ \\
426 & $0.280 \pm 0.019$  & $0.0073 \pm 0.0006 $ \\
600 & $0.376 \pm 0.026$  & $0.0097 \pm 0.0011 $ 
\end{tabular}
\end{ruledtabular}
\caption{\color{red}
\label{table:pt-energy-rms-and-drift} Standard deviation (SD) and drift of the total energy over 20~ps leapfrog trajectories, averaged over all trajectories at selected temperatures.}
\end{table}
Previous work has demonstrated that a Markov model based on a six-state decomposition, as depicted in Fig.~\ref{fig:energyLandscape}, can accurately describe the dynamics of this peptide for lagtimes longer than $\tau=6 \text{ ps}$~\cite{Chodera:2006p2352}. We employ the same state decomposition for all temperatures with the suggested minimal lagtime of $\tau=6 \text{ ps}$.  

\color{red}
To evaluate the accuracy of the methods for estimating transition probabilities, we compare the separate estimates obtained using dynamical reweighting (DR), transition matrix estimation (TE), and rate matrix estimation (RE) with a separate dataset of $6~\times~10~000$ short (10~ps) trajectories (ST) initiated from the equilibrium ensemble within each state at 302~K, also taken from Ref.~\cite{Chodera:2006p2352}. 
The PT simulation, in comparison, furnishes a total of 501 independent trajectories at that temperature.
\color{black}

The system is small enough that reasonable statistics can be obtained with moderate CPU requirements, while complex enough that some transitions (and even some states) are sampled only at high temperatures. 
\color{red}
In what follows, the results from the Markov model obtained from the dynamical reweighting method are compared to the model computed by Bayesian analyses using data from \emph{individual} temperatures from the parallel tempering simulation, as in Buchete and Hummer~\cite{Buchete:2008p2354}.
\color{black}

\begin{figure}[]
\includegraphics[width=8cm]{fig/plot_energyLandscapeAt302K.pdf}
\caption{\label{fig:energyLandscape} {\bf Terminally-blocked alanine peptide and potential of mean force with Markov state definitions.}
Left: The terminally blocked alanine peptide with $(\phi,\psi)$ torsions labeled. Right: The potential of mean force as a function of $(\phi,\psi)$ torsions at 302~K in units of $k_{B} T$, estimated from the PT simulation using WHAM~\cite{Chodera:2007p2346,Kumar:1992p860}, truncated at 10~$k_{B}T$ (white regions). The six manually identified states are labeled in black (note periodic boundaries)~\cite{Chodera:2006p2352}.
}
\end{figure}

\subsection{Estimated transition probabilities as a function of temperature}

A comparison of various approaches to estimating the transition probabilities between all $6\times 6$ pairs of states as a function of temperature is given in Fig.~\ref{fig:transition-probabilities-temp}. The blue lines give the estimates from dynamical reweighting (DR)~\cite{Chodera:2006p2352} using all available data at all temperatures, as described in Section~\ref{section:theory}.
To obtain the normalization constants $\hat{Z}_{k}$, we solved the set of self-consistent equations in Eq.~\ref{eq:free-energies} with a relative convergence tolerance in the residual of $10^{-7}$ (see Appendix~\ref{appendix:solving-self-consistent-equations}). Transition probabilities were also estimated at one intermediate temperature between each pair of simulated temperatures. 

The red lines in Fig.~\ref{fig:transition-probabilities-temp} show transition probabilities for the reversible single temperature estimation of transition matrices (TE)~\cite{Noe:2008p2459}. For each of the 40 temperatures the sampler was run to collect a total of 10~000 samples. For the sampling of reversible rate matrices (RE) the same amount of data was collected using the sampling proposed in Ref.~\cite{Buchete:2008p2354}. Diagnostics of convergence for both methods appear as supplementary Fig.~\ref{fig:convergence-of-bayesian-methods}. For reasons of clarity, the performance of the rate matrix estimation (RE) is only shown in the detailed comparison plots in Fig.~\ref{fig:detailed-probabilities-302K} discussed in the subsequent section. The black cross-hair in Fig.~\ref{fig:transition-probabilities-temp} refers to the reference values (ST) at 302~K estimated from the shooting trajectories. 

Qualitatively, both methods agree, especially for transitions among highly populated states (1,2,3,4). However, the dynamical reweighting estimate, which uses the combined data from all temperatures, has smaller uncertainties than the estimators that use only individual temperatures. In addition, the general agreement with the reference simulation is best for dynamical reweighting. 
%See Table~\ref{table:comparison-302K} for a detailed comparison at 302~K. % JDC: Table III is mentioned in the next section.

\begin{figure*}[!htb]
\includegraphics[width=\textwidth]{fig/plot_allTransitionProbabilites_woLimits.pdf}
\caption{\label{fig:transition-probabilities-temp} {\bf Comparison of all inter-state transition probabilities at a lagtime of $\tau=6 \text{ ps}$ as a function of temperature.}
Shaded regions denote $95\%$ confidence intervals about the estimated mean. 
Blue lines show the transition probabilities estimated using dynamical reweighting (DR). Red lines show the estimates from transition matrix estimation (TE) computed from only single temperature data. The black cross-hair indicates the reference using the shooting trajectory data (ST) at 302~K only.}
\end{figure*}

\subsection{Detailed comparison of transition probability estimates at 302 K}

For a detailed comparison with precisely known transition probabilities, the Bayesian analysis method with reversibility constraint for transition matrices~\cite{Noe:2008p2459} (TE) was also applied to a large set of shooting trajectories at 302~K, in which the trajectories were initiated from an equilibrium distribution within each state. \markRed{The results of this comparison at 302 K at a lagtime of $\tau = 6~\text{ps}$ between the distributions of transition probabilities $T_{ij}(\tau)$ of the different estimation methods are shown in Fig.~\ref{fig:detailed-probabilities-302K}.} All colors are consistent with Fig.~\ref{fig:transition-probabilities-temp} and Table~\ref{table:methods-colors}.

Overall the reweighting method performs very well compared to the single-temperature estimates. Even transition probabilities that are sampled very poorly at 302~K (such as transitions involving high free energy states 5 and 6) are in good agreement with the reference values at 302~K. Table~\ref{table:comparison-302K} shows the standard deviation in the absolute difference of the estimation methods compared to the reference simulation (ST) using a lagtime of $\tau = 6~\text{ps}$. The dynamical reweighting (DR) estimates have a smaller deviation than both Bayesian methods for both high and low free-energy states.

For transitions that are not sampled at certain temperature ranges, the maximum-likelihood estimates obtained with the present reweighting method are close to zero (see Fig.~\ref{fig:transition-probabilities-temp}). Generally, for transition probabilities close to zero or unity, the asymptotic normal distribution assumed in the statistical error estimate of dynamical reweighting is a poor approximation to this highly asymmetric distribution and therefore tends to overestimate the true uncertainties in these cases (see Fig.~\ref{fig:detailed-probabilities-302K}). 
%Thus, we suspect that reweighting will significantly overestimate the uncertainty in these cases since this is derived from the locally estimated Hessian of the probability distribution.

\color{red}
Comparing the two Bayesian methods, we find they behave almost identically for transitions among states where many transition counts are observed, but differ for transitions in which few transition events were observed (those involving states 5 and 6). 
Recall that these methods utilize the same likelihood functions but different parameterizations and priors; the influence of this difference is expected to be most prominent when statistics are poor, which is exactly as observed here.
\color{black}

\begin{table}
\begin{ruledtabular}
\begin{tabular}{lrcl}
 & \multicolumn{3}{c}{RMSE at 302 K} \\ \cline{2-4}
 & (DR) & (TE) & (RE)\\ \hline
Low free-energy states (1,2,3,4) & 0.007 & 0.019 & 0.020 \\
High free-energy states (5,6) & 0.140 & 0.223 & 0.162 \\ \hline
All Transitions & 0.079 & 0.126 & 0.092
\end{tabular}
\end{ruledtabular}
\caption{\label{table:comparison-302K} Root-mean squared error (RMSE), computed from the absolute difference of transition probabilities for $\tau = 6~\text{ps}$ to the reference simulation dataset (ST) at 302~K for the three methods of Markov model estimation and high and low free-energy subsets of transitions. 
%Dynamical reweighting (DR) shows the smallest RMSE compared to the reference simulation (ST).
}
\end{table}

\begin{figure*}[!htp]
\includegraphics[width=\textwidth]{fig/plot_allComparisonDetailed302K_woLimits.pdf}
\caption{\label{fig:detailed-probabilities-302K} \markRed{{\bf Detailed comparison of the distribution of transition probabilities $T_{ij}(\tau)$ for $\tau = 6~\text{ps}$ at 302~K for different estimation methods.} 
Red: Single temperature estimation of transition matrix (TE); 
Green: Single temperature estimation of rate matrix (RE); 
Blue: dynamical reweighting estimation (DR), Black: reference using shooting trajectories (ST). }}
\end{figure*}

\subsection{Comparison of temperature dependence of eigenvalues}
\markRed{
Dynamical reweighting can also be applied to estimate properties derived from the transition probabilities. For example, the eigenvalues $\lambda_{i}$ of a transition matrix are related to the timescales of processes $t_{i}^{*}$ by
\begin{equation}
	t_{i}^{*} = - \tau / \ln(\lambda_{i}),
	\label{eq:implied-timescale}
\end{equation}
where we assume that the eigenvalues $\lambda_{i}$ are sorted in order of descending absolute value ($ \lambda_{1} = 1 > | \lambda_{2} | > \ldots > |\lambda_{M}| > 0$)~\cite{Swope:2004p1280}. 
Note that Eq.~\ref{eq:implied-timescale} implies that eigenvalues close to unity are related to slow processes---those we are generally most interested in.
}

We investigated the dependence of the eigenvalues on the temperature in the present system. Fig.~\ref{fig:eigenvalues-temp} compares estimates for the second and third eigenvalues ($\lambda_{2}, \lambda_{3}$) of the transition matrix at each temperature with the different methods. The variance in the (TE) case was estimated from the set of eigenvalues of each sampled transition matrix. 
\color{red}
To estimate the statistical error in the estimates produce by (DR), we used a first-order Taylor expansion to propagate the statistical uncertainties in the transition matrix to uncertainties in the eigenvalues~\cite{Hinrichs:2007p5942}.
\color{black}
\begin{figure}[]
\includegraphics[width=8cm]{fig/plot_Eigenvalue_Vs_Temperature_No2.pdf}
\includegraphics[width=8cm]{fig/plot_Eigenvalue_Vs_Temperature_No3.pdf}
\caption{\label{fig:eigenvalues-temp} {\bf Temperature dependence of estimated eigenvalues.}
Red: Single temperature estimation of transition matrix (TE), Blue: dynamical reweighting estimation (DR), Black: reference 2nd and 3rd eigenvalue at 302~K using shooting trajectories (ST). Left: Comparison of the second largest eigenvalue vs temperature, Right: Comparison for third eigenvalue. The third reference eigenvalue is well predicted by both estimation methods at all temperatures although it matches only the second eigenvalue in the transition matrix estimation (TE). The second reference eigenvalue at low temperatures (below 350~K) is only detected by dynamical reweighting (DR).}
\end{figure}
At low temperatures (below 350~K), the second eigenvalue is estimated correctly by dynamical reweighting (DR), but not by single-temperature estimations. This is due to the fact that the transition process corresponding to this slowest timescale is not sampled at these low temperatures. Thus, estimates using only data collected at that temperature are erroneous. The agreement of dynamical reweighting timescales with the reference simulation is very good, although the error bars of the reweighted estimate are still very large compared to the good agreement of the estimated values with the reference values from the shooting trajectories. We speculate that the inappropriate approximation of the asymmetric posterior distributions with normal distributions used for the linear error propagation may lead here as well to an overestimation of the errors in the transition probabilities.

The third largest eigenvalue is predicted by both methods equally well (Fig.~\ref{fig:eigenvalues-temp}), although it occurs as the second-largest eigenvalue in the single-temperature estimates, which missed sampling the slowest process (described by $\lambda_{2}$) completely. A direct comparison of the predicted eigenvectors (Fig.~\ref{fig:eigenvectorMatchingComparison}) reveals that the slowest process (given by the second eigenvector of the reference transition matrix (ST)) is not detected by any of the single temperature methods. However, dynamical reweighting successfully finds all the processes, although the matching eigenvalues, and thus timescales, are permuted for faster processes. 

The comparison of Markov models is a nontrivial task~\cite{Bowman:2010p6510}, for which we use a symmetrized form of the transition matrix $\bfm{T}^{sym}$ and expand it into a sum of rank one matrices spanned by an outer product of the eigenvectors of $\bfm{T}^{sym}$ by
\begin{eqnarray}
	\bfm{T}^{sym} & = & \diag(\bfv{\pi}^{1/2}) \; \bfm{T} \; \diag(\bfv{\pi}^{-1/2}) \\
	              & = & \bfm{R} \; \diag(\bfv{\pi}^{1/2}) \; \bfm{R^{T}} \\
	              & = & \sum \limits_{i=1}^{M} \lambda_{i} \bfv{r}_{i} \; \bfv{r}_{i}^{T}.
\end{eqnarray}
\markRed{
Here, $\diag(\bfv{\pi}^{1/2})$ means the diagonal matrix with the square root of the equilibrium distribution $\bfv{\pi} = \{ \pi_{1}, \ldots, \pi_{M}\}$ on the diagonal, $\bfm{R} = \{ \bfv{r_{1}}, \ldots, \bfv{r_{M}} \}$ the matrix of normalized eigenvectors of $\bfm{T}^{sym}$ and $\lambda_{i}$ the corresponding eigenvalues. 
Neglecting the timescales (i.e.~the eigenvalues $\lambda_{i}$), two Markov models are similar if their inherent processes, described by the right eigenvectors $\bfv{r}_{i}$, are similar. Thus, we define the similarity matrix $\bfm{S}$ for two transition matrices $\bfm{T}$ and $\bfm{T}'$ by the mutual scalar product of their eigenvectors $\bfv{r}^{T}_{i}$ and $\bfv{r}'_{j}$ by
\begin{equation}
	S_{ij}\left[ \bfm{T}, \bfm{T}' \right] = \bfv{r}^{T}_{i} \bfv{r}'_{j}
\end{equation}
the results of which are presented in Fig.~\ref{fig:eigenvectorMatchingComparison}.
}
\begin{figure*}[!htb]
\includegraphics[width=\textwidth]{fig/plot_Eigenvector_Matching_Comparison.pdf}
\caption{\label{fig:eigenvectorMatchingComparison} \markRed{ {\bf Similarity matrices $\bfm{S}$ (scalar product) of eigenvectors from symmetrized transition matrices estimated with different methods at 302 K.} 
The eigenvectors indicate the states involved in the process, and thus high similarity (black) indicates a good approximation to the reference process (ST). Eigenvectors are sorted as descending eigenvalues. The 2nd eigenvector is found most reliably by dynamical reweighting (DR).}}
\end{figure*}
\subsection{Contributions from different temperatures to the estimates of expectation values}
\markRed{
The relative contribution from each temperature to the estimation of any expectation value at a given temperature is presented in Fig.~\ref{fig:contribution-to-counts}. 
We plot the relative total contribution from the subset of trajectories $\mathcal{Q}_{k}$ sampled from the distribution at $\beta_{k}$ to estimates at $\beta_{l}$ given by
\begin{equation}
	\bar{w}(\beta_{l} | \beta_{k}) = \frac{\sum \limits_{n \in \mathcal{Q}_{k}} w_{n}(\beta_{l})}{\sum \limits_{n \in \mathcal{Q}_{l}} w_{n}(\beta_{l})}
	\label{eg:average_weight}
\end{equation}
where the normalized trajectory weights $w_{n}(\beta)$ are defined in Eq.~\ref{equation:trajectory-weights}. On average, seven temperatures contribute more than 1\% to the expectation. 

Fig.~\ref{fig:contribution-to-weights} illustrates the contribution from the sampled data at $\beta_{k}$ to the computed expectation of the time-correlation function (Eq.~\ref{equation:definition-of-correlation}) for the rarely sampled transition $6 \rightarrow 5$.
This can be written as,
\begin{equation}
	\hat{C}_{65}(\beta_{l}\mid \beta_{k}) = \sum \limits_{n \in \mathcal{Q}_{k}} w_{n}(\beta_{l}) \hat{C}_{65}^{(n)},
\end{equation}
with $k$ indicating the temperature contributed from and $l$ the temperature estimated at. The weighted combination of estimations from multiple temperatures provides a much smoother and continuous estimation (rightmost column) then the estimation from a single temperature (topmost column). 
\begin{figure}[]
\includegraphics[width=7cm]{fig/plot_contribution_to_summary.pdf}
\caption{\label{fig:contribution-to-counts}{\bf Relative contribution $\bar{w}_{kl}$ in Eq.~\ref{eg:average_weight} to the estimates at inverse temperature $\beta_{k}$ from simulations at inverse temperature $\beta_{l}$ summed over all trajectories at the same temperature.} 
On average, seven temperatures contribute more than 1\% each to the estimation. }
\end{figure}
\begin{figure}[]
\includegraphics[width=7cm]{fig/plot_contribution_to_estimate_W65kl.pdf}
\caption{\label{fig:contribution-to-weights}{\bf Contributions to the estimation of the correlations $\hat{C}_{65}$ for the transition $6\rightarrow 5$ as defined in Eq.~\ref{equation:definition-of-correlation}.} 
The sum of one row (rightmost column), equal to the total counts estimated by the method at the desired temperature, provides a smoother estimate than the single temperature estimation (topmost column).}
\end{figure}
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DISCUSSION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Discussion}
\label{section:discussion}

\color{red}
We have presented a method that provides an estimate of Markov state transition probabilities from parallel tempering (PT) molecular dynamics simulations of biomolecules as a continuous function of temperature.
This allows Markov models to be constructed for intermediate temperatures not included in the simulation.
In addition, transition matrix estimates at temperatures included in the simulation are much more precise than those obtained with the methods examined here that make use of data from a single temperature alone. 
At low temperatures, even when some transitions are not observed at all at the temperature of interest, corresponding transition probabilities can still be estimated by incorporating dynamical information from higher temperatures without resorting to approximate rate laws such as Arrhenius. 
\color{black}

Additionally, the estimates of transition probabilities can be differentiated with respect to the inverse temperature $\beta$, because the trajectory weights $w_{n}(\beta)$ in Eq.~\ref{equation:trajectory-weights} are differentiable functions of temperature. 
This allows, in principle, kinetic and thermodynamic properties that depend on temperature derivatives (e.g., heat capacities) to be computed, provided care is taken in dealing with numerical issues since the trajectory weights $w_{n}(\beta)$ can easily span hundreds of orders of magnitude. 

\color{red}
In our illustrative calculations, we chose to employ a modified PT protocol to produce a series of Hamiltonian trajectories with initial configurations drawn from the NVT ensemble. 
%The requirement of good exchange probabilities among neigh ensures good overlap in the contribution to the dynamical reweighting for neighboring temperatures so that dynamical reweighting can be applied.  % JDC: This isn't exactly correct.  See my dynamical reweighting theory paper Discussion covering the generalized heat capacity.
However, the approach itself is not limited to Hamiltonian trajectories, but can be extended to other dynamics (e.g.~Brownian and Langevin dynamics) provided a temperature-independent dynamical analogue of the Hamiltonian can be computed, as described in Ref.~\cite{Chodera:2009p6073}.

The way in which the transition probabilities are estimated in terms of equilibrium correlation functions requires that the trajectory segments sampled during the PT simulation are drawn from equilibrium, which is ensured by the modified PT protocol provided the simulation is sufficiently long (see Appendix~\ref{appendix:parallel-tempering-stationary-distribution}). 
For systems with long correlation times, there may be very few \emph{uncorrelated} trajectories sampled from global equilibrium, but mixing within the Markov states may be sufficiently fast for many uncorrelated phase space configurations to be generated \emph{within} each state.
In this case, it is conceivably advantageous to apply dynamical reweighting to the trajectories originating from each set $\Gamma_{i} \subset \Omega$ separately to estimate each row of the transition matrix $\bfm{T}(\tau)$ separately, though the resulting matrix is no longer guaranteed to satisfy detailed balance.
\color{black}

%In addition, the validity of the estimation of the statistical error in all methods can be verified by using the  estimated statistical inefficiency to ensure  to reduce the data to an uncorrelated subset~\cite{Chodera:2009p6073}.

\markRed{The degree to which the use of PT can enhance the thermodynamic sampling efficiency is limited. Although activated processes will be sampled more often at higher temperatures, it becomes less that entropic bottlenecks will be penetrated. 
Lower temperatures, on the other hand, lead to an increased sampling of entropic barriers, while at the same time decreasing transitions across enthalpic barriers. 
This effectively limits the possible improvement in sampling and the range of temperatures that contribute substantial weight to transitions at a given temperature of interest. 
Despite this, dynamical reweighting allows information about activated processes is transferred from higher to lower temperatures and, for entropic barriers, from lower to higher temperatures. 
%In our illustrative example, transitions in the alanine peptide test system are dominated by activated processes (i.e., enthalpic barriers). % JDC: I'm not sure we can say that.
}

\markRed{
Both single-temperature methods give similar results for transitions with good statistics, differing mostly for transitions that were only rarely or not at all sampled, due presumably to the a dependence on the choice of Bayesian prior. The rate matrix estimation (RE) assumes, in addition to the detailed balance constraint, positivity of all transition matrix eigenvalues and non-negative off-diagonal rate matrix entries. The uniform distribution of parameters in logarithmic space leads most likely to favoring of low transition probabilities in states with poor transition statistics.
Surprisingly, the Bayesian estimates provide a reasonable bound on transition probabilities to and from a state even when the state is not sampled at all. This is most likely due to the reversibility constraint, which seems to provide information even in cases where there are few transitions to or from a state.
}

\markRed{
The general predictions of of dynamical reweighting are very good, while the quality of the statistical error estimation is limited near extremely small or large transition probabilities due to the reliance on asymptotic normality in the errors. 
Some combination of Bayesian and reweighting methods (such as T-WHAM~\cite{gallicchio:2005a}) may provide the best of both types of estimators by yielding more accurate uncertainties at the expense of introducing some bias from the introduction of energy histograms or some other parametric distributions for describing the energy density of states (now transitions). 
Finally, the enhanced estimates of mean values and their respective statistical uncertainties may be used to guide subsequent (potentially adaptive) sampling strategies, as described in Ref.~\cite{Singhal:2005p2464}.
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ACKNOWLEDGMENTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Acknowledgments}

The authors would like to thank Jed W. Pitera (IBM Almaden), Nicolae-Viorel Buchete (UCD Dublin), and Gerhard Hummer (NIH) for stimulating conversations during the execution of this work. 
JHP gratefully acknowledges funding from the German Research Foundation (DFG) through the award of a doctoral scholarship in the International Graduiertenkolleg IGK 710: ``Complex processes: Modeling, Simulation and Optimization''. 
JDC gratefully acknowledges support from HHMI and IBM predoctoral fellowship programs, NIH grant GM34993 through Ken A. Dill (UCSF), and NSF grant for Cyberinfrastructure (NSF CHE-0535616), and a California Institute for Quantitative Biosciences (QB3) Distinguished Postdoctoral Fellowship at various points throughout this work.
VSP acknowledges support from NIH RO1 GM062868.
JHP and FN acknowledge support from DFG Research Center Matheon and DFG grant No. 725/2. JCS acknowledges funding from the U.S. Department of Energy ``Multiscale Mathematics'' SciDAC and Genomes-to life program under grant No. ERKJE84/ERKPE84.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BIBLIOGRAPHY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\bibliographystyle{prsty} 

\begin{thebibliography}{10}

\bibitem{kern:science:2002}
E.~Z. Eisenmesser, D.~A. Bosco, M. Akke, and D. Kern, Science {\bf 295},  1520
  (2002).

\bibitem{kern:nature:2005}
E.~Z. Eisenmesser, O. Millet, W. Labeikovsky, D.~M. Korzhnev, M. Wolf-Watz,
  D.~A. Bosco, J.~J. Skalicky, L.~E. Kay, and D. Kern, Nature {\bf 438},  117
  (2005).

\bibitem{Feng:2005p6514}
H. Feng, Z. Zhou, and Y. Bai, Proc. Nat. Acad. Sci. USA {\bf 102},  5026
  (2005).

\bibitem{Min:2005p4264}
W. Min, G. Luo, B. Cherayil, S. Kou, and X. Xie, Phys. Rev. Lett {\bf 94},  1
  (2005).

\bibitem{Smith:2008p4287}
G. Smith, K. Lee, X. Qu, Z. Xie, J. Pesic, T. Sosnick, T. Pan, and N. Scherer,
  J. Mol. Biol. {\bf 378},  941  (2008).

\bibitem{zhuang:2003a}
X. Zhuang and M. Rief, Curr. Opin. Struct. Biol. {\bf 13},  88  (2003).

\bibitem{Matagne:1997p6517}
A. Matagne, S. Radford, and C.~M. Dobson, J. Mol. Biol. {\bf 267},  1068
  (1997).

\bibitem{Neuweiler:2005p5935}
H. Neuweiler, S. Doose, and M. Sauer, Proc. Nat. Acad. Sci. USA {\bf 102},
  16650  (2005).

\bibitem{Goldbeck:1999p6516}
R. Goldbeck, Y. Thomas, E. Chen, R. Esquerra, and D. Kliger, Proc. Nat. Acad.
  Sci. USA {\bf 96},  2782  (1999).

\bibitem{geyer-thompson:j-am-stat-assoc:1995:parallel-tempering}
C.~J. Geyer and E.~A. Thompson, J. Am. Stat. Assoc. {\bf 90},  909  (1995).

\bibitem{hukushima-nemoto:j-phys-soc-jpn:1996:parallel-tempering}
K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. {\bf 65},  1604  (1996).

\bibitem{hansmann:chem-phys-lett:1997:parallel-tempering-monte-carlo}
U.~H.~E. Hansmann, Chem. Phys. Lett. {\bf 281},  140  (1997).

\bibitem{Sugita_CPL314_141}
Y. Sugita and Y. Okamoto, Chem. Phys. Lett. {\bf 314},  141  (1999).

\bibitem{Jager:2001p6248}
M. J{\"a}ger, H. Nguyen, J.~C. Crane, J.~W. Kelly, and M. Gruebele, J. Mol.
  Biol. {\bf 311},  373  (2001).

\bibitem{Buchete:2008p2351}
N.-V. Buchete and G. Hummer, J. Phys. Chem. B {\bf 112},  6057  (2008).

\bibitem{Noe:2007p867}
F. No{\'e}, I. Horenko, C. Sch{\"u}tte, and J.~C. Smith, J. Chem. Phys. {\bf
  126},  155102  (2007).

\bibitem{Chodera:2006p2352}
J.~D. Chodera, W.~C. Swope, J.~W. Pitera, and K.~A. Dill, Multiscale Model.
  Sim. {\bf 5},  1214  (2006).

\bibitem{Muff:2009p6071}
S. Muff and A. Caflisch, J. Phys. Chem. B {\bf 113},  3218  (2009).

\bibitem{Chodera:2007p2463}
J.~D. Chodera, N. Singhal, V.~S. Pande, K.~A. Dill, and W.~C. Swope, J. Chem.
  Phys. {\bf 126},  155101  (2007).

\bibitem{Hinrichs:2007p5942}
N.~S. Hinrichs and V.~S. Pande, J. Chem. Phys. {\bf 126},  244101  (2007).

\bibitem{Singhal:2005p2464}
N. Singhal and V.~S. Pande, J. Chem. Phys. {\bf 123},  204909  (2005).

\bibitem{Swope:2004p907}
W.~C. Swope, J.~W. Pitera, and F. Suits, J. Phys. Chem. B {\bf 108},  6571
  (2004).

\bibitem{Swope:2004p1280}
W.~C. Swope, J.~W. Pitera, F. Suits, M. Pitman, M. Eleftheriou, B.~G. Fitch,
  R.~S. Germain, A. Rayshubskiy, T.~C. Ward, and Y. Zhestkov, J. Phys. Chem. B
  {\bf 108},  6582  (2004).

\bibitem{sarich:errorindiscretization}
M. Sarich, F. No{\'e}, and C. Sch{\"u}tte, Multiscale Model. Sim.  (2009).

\bibitem{Buchete:2008p2354}
N.-V. Buchete and G. Hummer, Phys. Rev. E {\bf 77},  4  (2008).

\bibitem{Ferrenberg:1989p3671}
A.~M. Ferrenberg and R.~H. Swendsen, Phys. Rev. Lett. {\bf 63},  1195  (1989).

\bibitem{bartels:1997a}
C. Bartels and M. Karplus, J. Comput. Chem. {\bf 18},  1450  (1997).

\bibitem{Kumar:1992p860}
S. Kumar, D. Bouzida, R. Swendsen, P. Kollman, and J. Rosenberg, J. Comp. Chem.
  {\bf 13},  1011  (1992).

\bibitem{gallicchio:2005a}
E. Gallicchio, M. Andrec, A.~K. Felts, and R.~M. Levy, J. Phys. Chem. B {\bf
  109},  6722  (2005).

\bibitem{Bartels:2000p4835}
C. Bartels, Chem. Phys. Lett. {\bf 331},  446  (2000).

\bibitem{Shirts:2008p2350}
M.~R. Shirts and J.~D. Chodera, J. Chem. Phys. {\bf 129},  124105  (2008).

\bibitem{Chodera:2009p6073}
J.~D. Chodera, M.~R. Shirts, J.-H. Prinz, W.~C. Swope, F. No{\'e}, and V.~S.
  Pande, J. Chem. Phys., accepted (2010).

\bibitem{Mitsutake:2001p3675}
A. Mitsutake, Y. Sugita, and Y. Okamoto, Biopolymers {\bf 60},  96  (2001).

\bibitem{Voter:1985p2344}
A.~F. Voter and J.~D. Doll, J. Chem. Phys. {\bf 82},  80  (1985).

\bibitem{Adams:1981p2340}
J. Adams and J.~D. Doll, Surf. Sci. {\bf 111},  492  (1981).

\bibitem{Chodera:2006p6224}
J.~D. Chodera, W. Swope, J. Pitera, and K.~A. Dill, Multiscale Model. Sim. {\bf
  5},  1214  (2006).

\bibitem{Earl:2005p2409}
D.~J. Earl and M.~W. Deem, Phys. Chem. Chem. Phys. {\bf 7},  3910  (2005).

\bibitem{andersen:1980:jcp:andersen-thermostat}
H.~C. Andersen, J. Chem. Phys. {\bf 72},  2384  (1980).

\bibitem{Noe:2008p2459}
F. No{\'e}, J. Chem. Phys. {\bf 128},  244103  (2008).

\bibitem{verlet:1967:phys-rev:verlet-integrator-1}
L. Verlet, Phys. Rev. {\bf 159},  98  (1967).

\bibitem{verlet:1967:phys-rev:verlet-integrator-2}
L. Verlet, Phys. Rev. {\bf 165},  201  (1968).

\bibitem{berendsen:1977:j-comput-phys:SHAKE}
J. Ryckaert, G. Ciccotti, and H.~J.~C. Berendsen, J. Comput. Phys. {\bf 23},
  327  (1977).

\bibitem{Chodera:2007p2346}
J.~D. Chodera, W.~C. Swope, J.~W. Pitera, C. Seok, and K.~A. Dill, J. Chem.
  Theo. Comp. {\bf 3},  26  (2007).

\bibitem{Bowman:2010p6510}
G.~R. Bowman, D.~L. Ensign, and V.~S. Pande, J. Chem. Theo. Comp. {\bf 6},  787
   (2010).

\bibitem{kinchin:statmech:1949}
A.~I. Kinchin, {\em Mathematical Foundations of Statistical Mechanics} (Dover
  Publications, NY, USA, 1949).

\end{thebibliography}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% APPENDICES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\appendix
\begin{widetext}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% APPENDIX: SOLVING THE SELF-CONSISTENT EQUATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{\label{appendix:solving-self-consistent-equations}
Efficient solution of the self-consistent equations for canonical distribution of Hamiltonian trajectories}

For the case of a canonical distribution of Hamiltonian trajectories, the normalization constants $\hat{Z}_{k}$ or alternatively the dimensionless free energies $\hat{f}_{i} \equiv - \ln \hat{Z}_{i}$ are defined through a set of $K$ coupled nonlinear equations
\begin{eqnarray}
\hat{f}_i &=& - \ln \sum_{n=1}^{N} \left[ \sum\limits_{k=1}^K N_{k} \, \exp[\hat{f}_k - (\beta_k - \beta_i) E_{n}] \right]^{-1} 
\label{equation:self-consistent-equation}
\end{eqnarray}
where all symbols as defined as for Eq.~\ref{eq:free-energies}. Any numerically stable method for solving a set of coupled nonlinear equations can, in principle, be used to obtain the $\hat{f}_i$.
A scheme for solving a more general form of these equations by self-consistent iteration or Newton-Raphson is described in Appendix C of~\cite{Shirts:2008p2350}.

Because of the structure of this specific case, we can rapidly obtain a close initial guess for the $\hat{f}_i$ by using a form inspired by the weighted histogram analysis method (WHAM)~\cite{Kumar:1992p860}.
By instead constructing $M$ bins in the total energy $E$ spanning a range $(E_\mathrm{min}, E_\mathrm{max})$, we can approximate Eq.~\ref{equation:self-consistent-equation} with a sum over histograms (as in Eqs. 19--20 of~\cite{Kumar:1992p860}):
\begin{eqnarray}
	\hat{f}^{(n+1)}_i &=& - \ln \sum_{m=1}^M H_m \, \left[ \sum_{k=1}^K N_k \, \exp[\hat{f}^{(n)}_k - (\beta_k - \beta_i) E_m \right]^{-1} 
\end{eqnarray}
where $H_m$ denotes the number of samples $E_{n}$ falling in histogram bin $m$, and $E_m$ represents the energy at the midpoint of that bin. For the number of bins, typically, a value of $M \approx 100$ can be used.
Since Eq.~\ref{equation:self-consistent-equation} is linear in the $\hat{Z}_{k}$, the $\hat{f}_{k}$ are unique up to an additive constant and we can fix one value, say $f_1$, by subtracting off the computed value of $f^{(n+1)}_1$ after each iteration in order to avoid numerical drift.

After the initial guess has been reached, self-consistent iteration can rapidly refine the free energies to the desired tolerance while eliminating the bias arising from the use of histograms:
\begin{eqnarray}
\hat{f}^{(n+1)}_i &=& - \ln \sum_{n=1}^{N} \left[ \sum\limits_{k=1}^K N_{k} \, \exp[\hat{f}^{(n)}_k - (\beta_k - \beta_i) E_{n}] \right]^{-1} 
\end{eqnarray}
Again, we fix $\hat{f}_1 = 0$ and terminate iterations when a relative tolerance $\max_{i = 2,\ldots,K} |\hat{f}^{(n+1)}_i - \hat{f}^{(n)}_i| / |\hat{f}^{(n+1)}_i + \hat{f}^{(n)}_i|$ is less than some given tolerance that ensures the computed expectations of properties of interest are no longer changing.
We find that $10^{-7}$ is often a safe choice.

Cautions observed in Appendix C of~\cite{Shirts:2008p2350} regarding sums of logarithms and numerical over/underflow in the evaluation of exponentials should be observed in implementation of this, or any, algorithm for obtaining the $\hat{f}_i$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% APPENDIX: PROOF THAT MODIFIED PARALLEL TEMPERING PROTOCOL GENERATES CORRECT DISTRIBUTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{\label{appendix:parallel-tempering-stationary-distribution}Proof that modified PT protocol generates canonical distribution}

Here, we prove that the modified PT protocol described in Section~\ref{section:modified-parallel-tempering-protocol} samples from the canonical stationary distribution at all temperatures.

Define stationary distributions for momenta $\p$ and coordinates $\q$ in Cartesian space $\mathbb{R}^{3N}$ at inverse temperature $\beta$:
\begin{eqnarray}
\pi_p(\p | \beta) = [P(\beta)]^{-1} \, e^{-\beta T(\p)} &;& P(\beta) = \int d\p \, e^{-\beta T(\p)} \nonumber \\
\pi_q(\q | \beta) = [Q(\beta)]^{-1} \, e^{-\beta U(\q)} &;& Q(\beta) = \int d\q \, e^{-\beta U(\q)} 
\end{eqnarray}
where $T(\p)$ denotes the kinetic energy and $U(\q)$ the potential energy function.
Suppose we have two replicas whose current phase space points are denoted by $\z_1 = (\q_1, \p_1)$ and $\z_2 = (\q_2, \p_2)$, initially at equilibrium at their respective inverse temperatures $\beta_1$ and $\beta_2$, such that
\begin{eqnarray}
\p_1 \sim \pi_p(\p_1 | \beta_1) &;& \q_1 \sim \pi_q(\q_1 | \beta_1) \nonumber \\
\p_2 \sim \pi_p(\p_2 | \beta_2) &;& \q_2 \sim \pi_q(\q_2 | \beta_2) .
\end{eqnarray}
We now consider what happens to the distributions of $\z_1$ and $\z_2$ after an exchange attempt.
Define ``post-exchange attempt'' coordinates and momenta for inverse temperature $\beta_1$:
\begin{eqnarray}
\q'_1 &\leftarrow& \begin{cases}
\q_1 \quad \mathrm{with\:prob.} \quad 1 - \theta(\q_1, \q_2 | \beta_1, \beta_2) & \mathrm{(rejected)} \\
\q_2 \quad \mathrm{with\:prob.} \quad \theta(\q_1, \q_2 | \beta_1, \beta_2) & \mathrm{(accepted)} \\
\end{cases} \nonumber \\
\p'_1 &\sim& \pi_p(\p'_1 | \beta_1) \quad \mathrm{(velocity\:randomization)} \nonumber
\end{eqnarray}
where the exchange acceptance probability $\theta(\q_1, \q_2 | \beta_1, \beta_2)$ is given by
\begin{eqnarray}
\theta(\q_1, \q_2 | \beta_1, \beta_2) &=& \min\{ 1, \exp[-\beta_1 U(\q_2) - \beta_2 U(\q_1) + \beta_1 U(\q_1) + \beta_2 U(\q_2) \}
\end{eqnarray}
We now compute the distribution of $\q'_1$, the configuration supposedly at temperature $\beta_1$ after the exchange attempt:
\begin{eqnarray}
\rho_1(\q'_1) &=& \int d\q_2 \, [1 - \theta(\q'_1, \q_2 | \beta_1, \beta_2)] \, \pi_q(\q'_1 | \beta_1) \, \pi_q(\q_2 | \beta_2) + \int d\q_2 \, \theta(\q_2, \q'_1 | \beta_1, \beta_2) \, \pi_q(\q_2 | \beta_1) \, \pi_q(\q'_1 | \beta_2) \nonumber \\
&=& \int d\q_2 \, [1 - \min\{ 1, e^{-\beta_1 U(\q_2)} \, e^{- \beta_2 U(\q'_1)} \, e^{+ \beta_1 U(\q'_1)} \, e^{+ \beta_2 U(\q_2)} \}] \, \frac{e^{-\beta_1 U(\q'_1)}}{Q(\beta_1)} \, \frac{e^{-\beta_2 U(\q_2)}}{Q(\beta_2)} \nonumber \\
& &\qquad + \int d\q_2 \, \min\{ 1, e^{-\beta_1 U(\q'_1)} \, e^{- \beta_2 U(\q_2)} \, e^{+ \beta_1 U(\q_2)} \, e^{+ \beta_2 U(\q'_1)} \} \, \frac{e^{-\beta_1 U(\q_2)}}{Q(\beta_1)} \, \frac{e^{-\beta_2 U(\q'_1)}}{Q(\beta_2)} \nonumber \\
&=& \frac{e^{-\beta_1 U(\q'_1)}}{Q(\beta_1)} - \int d\q_2 \, \min\left\{ \frac{e^{-\beta_1 U(\q'_1)}}{Q(\beta_1)} \, \frac{e^{-\beta_2 U(\q_2)}}{Q(\beta_2)} , \frac{e^{-\beta_1 U(\q_2)}}{Q(\beta_1)} \, \frac{e^{-\beta_2 U(\q'_1)}}{Q(\beta_2)} \right\} \nonumber \\
& &\qquad + \int d\q_2 \, \min\left\{ \frac{e^{-\beta_1 U(\q_2)}}{Q(\beta_1)} \, \frac{e^{-\beta_2 U(\q'_1)}}{Q(\beta_2)}, \frac{e^{-\beta_1 U(\q'_1)}}{Q(\beta_1)} \, \frac{e^{-\beta_2 U(\q_2)}}{Q(\beta_2)} \right\}  \nonumber \\
&=& \pi_q(\q'_1 | \beta_1)
\end{eqnarray}
Therefore, after the exchange attempt, the new configuration $\q'_1$ is still at equilibrium at inverse temperature $\beta_1$.
(A similar series of steps can be applied for the temperature $\beta_2$.)

Redrawing the momentum from the Maxwell-Boltzmann distribution at inverse temperature $\beta_1$ will, of course, not change the equilibrium distribution, and can be shown to only support the canonical distribution at inverse temperature $\beta_1$, and no other stationary distribution~\cite{andersen:1980:jcp:andersen-thermostat}.
Evolution by Hamiltonian dynamics for any length of time does not alter the stationary canonical distribution~\cite{kinchin:statmech:1949}.
Therefore, the proposed protocol samples from the canonical distribution at the desired temperatures, provided sufficient time is allowed for equilibration.  \qed

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% APPENDIX: PROPAGATION OF ERROR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Convergence of transition probabilities in Bayesian Methods}
The convergence of transition probabilities from the Bayesian sampling methods is presented in Supplementary Fig.~\ref{fig:convergence-of-bayesian-methods} for various temperatures. Therefore, the proposed protocol generates the canonical distribution at the desired temperature.
\begin{figure*}[]
\includegraphics[width=15cm]{fig/plot_ConvergenceWithSamplesST.pdf}
\\ \vspace{0.5cm}
\includegraphics[width=15cm]{fig/plot_ConvergenceWithSamplesHT.pdf}
\caption{\label{fig:convergence-of-bayesian-methods} 95\% confidence intervals of Transition Probabilities sampled by the transition matrix estimation (upper Plot) (TE) and rate matrix estimation (lower plot) (RE) versus number of drawn samples. Color indicates performance by temperature. Blue: 302~K, Yellow: 426~K, Red: 600~K. After about 5~000 samples the confidence intervals stabilize suggesting reasonably well sampled transition probabilities.}
\end{figure*}
\clearpage
\end{widetext}
\end{document}