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

\usepackage{palatino}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}
%\usepackage[mathbf,mathcal]{euler}
%\usepackage{citesort}
\usepackage{dcolumn}
\usepackage{boxedminipage}
\usepackage{verbatim}
\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{\bfm}[1]{{\mbox{\boldmath{$#1$}}}}
%\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
%\newcommand{\x}{\bfv{x}}
%\newcommand{\y}{\bfv{y}}
%\newcommand{\f}{\bfv{f}}
%\newcommand{\z}{\bfv{z}}
%\newcommand{\p}{\bfv{p}}
%\newcommand{\q}{\bfv{q}}
%\newcommand{\T}{\mathrm{T}}

\newcommand{\op}[1]{\mathcal{#1}}
%\newcommand{\vektor}[1]{{\bf #1}}
%\renewcommand{\vec}[1]{\bfm{#1}} % italicized boldface
%\renewcommand{\vec}[1]{{\bf #1}}
\renewcommand{\vec}[1]{{\mbox{\boldmath{$#1$}}}}
\newcommand{\q}{\vec{q}} % coordinates
\newcommand{\p}{\vec{p}} % momenta
\newcommand{\z}{\vec{z}} % phase space point
\newcommand{\x}{\vec{x}} % general vector x
\newcommand{\M}{\bfm{M}} % diagonal mass matrix
\newcommand{\grad}{\nabla}
\newcommand{\timeavg}[1]{\overline{#1}}                    % time average over a trajectory

\newcommand{\bfc}{\bfm{c}}
\newcommand{\hatf}{\hat{f}}

\newcommand{\bTheta}{\bfm{\Theta}}
\newcommand{\btheta}{\bfm{\theta}}
\newcommand{\bhatf}{\bfm{\hat{f}}}
\newcommand{\Cov}[1] {\mathrm{cov}\left( #1 \right)}
\newcommand{\Ept}[1] {{\mathrm E}\left[ #1 \right]}
\newcommand{\Eptk}[2] {{\mathrm E}\left[ #2 \,|\, #1\right]}
\newcommand{\T}{\mathrm{T}}                                % T used in matrix transpose
\newcommand{\qed}{\hfill \mbox{\raggedright \rule{.07in}{.1in}}}

\newcommand{\DeltaT}{\Delta t}
%\newcommand{\tauD}{\tau_{\Delta}}
\newcommand{\tauD}{\hat{\tau}}
%\newcommand{\LD}{L_{\Delta}}
\newcommand{\LD}{\hat{L}}

\DeclareMathOperator{\diag}{diag}

%\linespread{2.0}
%\singlespacing
%\onehalfspacing

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

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

\author{Jan-Hendrik Prinz}
 \email{jan-hendrik.prinz@iwr.uni-heidelberg.de}
% \thanks{Corresponding author (both authors have contributed equally to this paper)}
 \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}
 \thanks{Corresponding author}
 \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 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 parallel tempering 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 rates 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 dipeptide.
\end{abstract}

%\pacs{02.30.Cj,05.70.Ce,82.20.Wt}

\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 of central importance in biology and medicine. {\color{red}[WHY?]} 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}, fluorescence emission \cite{Min:2005p4264,Smith:2008p4287}, energy transfer\cite{zhuang:2003a}, correlation spectroscopy \cite{Neuweiler:2005p5935}, and non-equilibrium perturbation experiments\cite{Smith:2008p4287}. 
%[Lewis Kay's work on ``seeing the invisible'', Dorothee Kern];
%[Sunney Xie anti-fluorescein]
%[Zhaowei Zhuang??]
%[Norbert Scherer Mg2+ cycling, T-jump experiments with multiple timescales]
Developing a quantitative understanding of what gives rise to these conformational states and the interactions that govern their transitions will provide insight into how post-translational
%{\color{red} [JCS] Why post-translational?} 
modification and non-covalent association can affect dynamics and function, and will have a significant impact on our understanding of many biological processes, such as signaling events, enzyme regulation, allostery, and drug design for conformationally flexible molecules.

Integrating the dynamic in biomolecular systems to sample the underlying phase-space often suffers from the problem that the timescales of conformational changes are usually orders of magnitude larger than simulationtimes the available computational ressources allow for. In these cases parallel tempering molecular dynamics simulations have been an effective and thus popular approach to overcoming the issue of convergence in molecular simluations by allowing the bath temperature to change as the simulation proceeds\cite{Sugita_CPL314_141}. At the same time this approach permits an analysis of temperature dependence, which is especially interesting for comparisons with certain experimental results\cite{Jager:2001p6248} (e.g. melting curves, heat capacities). Although parallel tempering has the drawback that unphysical replica trajectories are produced, the short physical Newtonian trajectories in between these exchanges can provide useful dynamical information. If the parallel tempering simulation is well-equilibrated, these short trajectoriy segments will be sampled from equilibrium of their coresponding temperatures.

Markov models provide a way to model the slow conformational dynamics of biomolecules such as peptides and proteins based on short simulations \cite{Buchete:2008p2351,Noe:2007p867,Chodera:2006p2352,Muff:2009p6071}.
%{\color{red} [JHP:] More Citations?}
In these models, generally, conformational states are envisioned as disjoint but connected regions of configurational space. The biomolecule spends long times within individual regions before undergoing rapid stochastic transitions between regions. 
If a separation of timescales exists between fast motions within a region and due to slow motions (i.e. infrequent transitions) between regions, the inter-state dynamics can be well described by a Markov Model in discrete timesteps $\tau$, where the coarse graining in time is necessary since the discretization in space prohibits us from trading arbitrarily fast relaxation processes. 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 shown recently that $\tau$ can be reduced (and the quality of the model is enhanced) by using additional states which are no longer metastable\cite{sarich:errorindiscretization}.

%It has been shown using appropriate definitions of the discrete states that the essential dynamics of biomolecules can be well-described by such a Markov model\cite{Buchete:2008p2351,Chodera:2006p2352,Noe:2007p867}. Still the construction of this model requires several components: 

%First, regions of conformation space that are metastable with respect to some timescale of interest must be enumerated. [cite state decomposition and spectral clustering papers and theses] Next, the transition probabilities or rates between states must be estimated. This is generally done by computing time-correlation functions \cite{Swope:2004p907} or counting [cite]. Finally, the model has to be validated: It must be verified that the model can recapitulate the observed kinetics to, at minimum, show that it is self consistent. [cite]

Recently, Buchete and Hummer have shown that thermodynamics and kinetics can be modeled over the range of temperatures by constructing Markov models using the short physical trajectories generated from parallel tempering simulations\cite{Buchete:2008p2354}.
However, if a complete description of dynamics across the entire thermally-accessible configurational space is desired, one quickly runs into problems if this 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 temperature of interest\cite{Buchete:2008p2354}. 

Furthermore, 
%independently of issues of temperature dependence, the configurational space spanned over the full range of temperatures simulated is larger than that at any single temperature. Hence, 
one would like to make use of the data collected at all temperatures to characterize the kinetic behavior in all regions sampled over the full range of temperatures spanned by the parallel tempering simulation. This approach resembles that of equilibrium reweighting.

Here, we propose an approach for integrating data from all temperatures by making use of dynamical reweighting\cite{Shirts:2008p2350}, allowing a continuous 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 multiple histogram reweighting\cite{Ferrenberg:1989p3671,Nose:1984p5958}, WHAM\cite{Kumar:1992p860,gallicchio:2005a} and MBAR\cite{Shirts:2008p2350}) allow the use of samples collected from multiple distributions to provide an improved estimate of the expectation value of some 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}.
%{\color{red} [JHP:][cite Sugita Okamoto papers. Not sure if this is a good one or we need more!}.

Dynamical reweighting has recently been proposed as a way of estimating dynamical properties (such as correlation functions) using an asymptotically optimal estimator, and also provides an 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) as a smooth function of temperature, making use of data from all temperatures. This has the advantage of producing estimates of transition probabilities at \emph{any} temperature and permitting the computation of properties that depend on derivatives with respect to temperature (e.g. heat capacities), along with an assessment of the error in the quantities.

We illustrate this approach for the standard test case of the terminally-blocked alanine dipeptide 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}. The dipeptide system presents a challenge for estimators based on individual temperatures, because states exist with very high free energies relative to the most populated states. These are poorly sampled at temperatures near 300 K, even though they dominate the long-time 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 parallel tempering simulation.

%\subsubsection{Outline of this manuscript}

This paper is organized as follows:
First we review the theory behind Markov models of multistate conformational dynamics. We then show how dynamical reweighting can be used to estimate temperature-dependent transition probabilities and rates for a given state decomposition. Finally, we illustrate the method by applying it to a six-state decomposition of an MD simulation of the terminally-blocked alanine dipeptide, and compare the results to the approach of Buchete and Hummer, in which Bayesian estimates of transition probabilities are obtained from a single temperature alone.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THEORY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theory}
\label{section:theory}
\subsection{Markov Models}

Let $\Omega$ be the configuration space with a complete decomposition into $M$ disjoint sets $\Gamma_{i} \in \Omega$ such that
\begin{equation}
	\bigcup\limits_{i=1}^{M} \Gamma_{i} = \Omega, \; \Gamma_{i} \cap \Gamma_{j} = \emptyset, \; \forall i, j \in \{1,\ldots,M\}, \; i \neq j.
\end{equation}
For convenience we also define indicator functions $\chi_{i}(q\bfv{q}) \in \{0,1\} $ by setting
\begin{equation}
	\chi_{i}(\bfv{q}) = 
	\begin{cases}
		1 & \text{if} \; \bfv{q} \in \Gamma_{i} \\
		0 & \text{else}
	\end{cases}
\end{equation}
i.e., the function assumes the value of unity if $\bfv{q}$ belongs to set $i$, and zero otherwise. Based on this discretization of state space, we can define a row-stochastic transition matrix $\bfm{T}$ with conditional probabilities of finding the system in state $j$ t time $\tau$ after it was originally in state $i$:
\begin{eqnarray}
T_{ij}(\tau) & = & \mathbb{P}\left(\bfv{q}(t+\tau) \in \chi_{j} \mid \bfv{q}(\tau)\right) \\
 & = & \frac{\left\langle \chi_{i}(\bfv{q}(t)) \chi_{j}(\bfv{q}(t+\tau)) \right\rangle}{\left\langle \chi_{i} \right\rangle}, \label{equation:transition-matrix-definition}
\end{eqnarray}
%\begin{eqnarray}
%	T_{ij}(\tau) &:=& \mathbb{P}(q(t+\tau) \in \Gamma_{j} | q(t)  \in \Gamma_{i}). \label{equation:transition-matrix-definition}
%\end{eqnarray}
which then evolves a probability distribution $\bfv{p} \in [0,1]^{M}$ by time step $\tau$ by 
\begin{equation}
	\bfv{p}(t+\tau) = \bfv{p}(t) \bfm{T}(\tau).
\end{equation}
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}
Eq. \ref{eq:approximateProbabilityTransport} is only an approximation due to the introduction of the coarse-graining $\Gamma \equiv \left\{ \Gamma_{1}, \ldots, \Gamma_{M}\right\}$. It has been shown \cite{sarich:errorindiscretization}, that the approximation error by the discretization $\Gamma$ can be made arbitrarily small by either choosing more states $M$ or increasing the underlying lagtime $\tau$.
%, provided, that the underlying real dynamics is Markovian in the full continuous space. 
Thus, with an appropriate choice of states we can ensure an approximation of the true dynamics by a time- and space-discrete Markov model to the necessary precision. Although this is a crucial step, the process of finding an optimal decomposition of state space is beyond the scope of this paper and has been discussed at length elsewhere \cite{Noe:2007p867,Chodera:2007p2463}.

\subsection{Estimating transition probabilities vs estimating transition rates}
In many studies coarse-grained dynamics is equivalently described by the \emph{master equation}
\begin{equation}
	\dot{\bfv{p}}(t) = \bfv{p}(t) \; \bfm{K}
\end{equation}
where $\bfm{K} \in \mathbb{R}^{m\times m}$ denotes the rate matrix \cite{Buchete:2008p2354,Voter:1985p2344}{\color{red}[CITE: Straub & Doll ?]}. While $\bfm{T}(\tau)$ can be straightforwardly estimated from a series of observed transitions at time step $\tau$ using Eq. \ref{equation:transition-matrix-definition}, $\bfm{K}$ cannot because the inversion of the equation
\begin{equation}
	\bfm{T}(\tau) = \exp( \tau \bfm{K} )
\end{equation}
is not unique for stochastic matrices unless $\bfm{T}(\tau)$ is positive definite and reversible, but even then represents the correct dynamics only for times larger than some Markov time $\tau^{*}. In the following we focus on the estimation of time-discrete transition probabilities{\color{red} CITE: Hummer papers about rate matrix inference}.
%While the optimal estimate of transition probabilities from dynamical reweighting gives a transition matrix $\bfm{T}(\tau)$ that satisfies detailed balance, it is nontrivial to compute a valid rate matrix $\bfm{K}$ from $\bfm{T}$. In theory, $\bfm{K}$ and $\bfm{T}$ are related by \ref{eq:ratematrix-to-transitionmatrix}, but in practice, statistical and numerical errors in the estimated matrix $\bfm{T}$ that produce negative (though generally small) eigenvalues, frustrating the calculation of the matrix logarithm. Even when $\bfm{T}$ possesses no negative eigenvalues, additional constraints must be satisfied to produce a $\bfm{K}$ with positive off-diagonal rates. While an optimization scheme could find a suitable matrix $\bfm{K}$ such that the error $||\bfm{T}(\tau) - \exp[\bfm{K}]||$ is minimized, we consider the details of such a process and the optimal way to do so beyond the scope of this paper.

\subsection{Estimation from time-discrete trajectories}
%{\color{red}[JHP]: Keep theory only time-continuous and mention, that it can be done time-discrete.}
%Obviously, our model must be constructed from finite statistics, which leaves us with the choice of an appropriate transition matrix to represent the observed dynamics. The theory can be written completely for time-continuous trajectories, but since our observations will be time-discrete we continue from here with discrete time-series of states as observations.

Consider a trajectory $\bfv{q}(t)$ sampled at time intervals $\DeltaT$ and let $\tauD \in \mathbb{N}$ be the lagtime $\tau = \tauD \Delta t$ and respectively $\LD \in \mathbb{N}$ the length of the trajectory in units of $the sampling interval \DeltaT$. We introduce $\tauD = \tau / \DeltaT \in \mathbb{N}$ and $\LD = L / \DeltaT \in \mathbb{N}$ and define a time discrete trajectory $\bfm{q} \in \Omega^{\LD+1}$
\begin{equation}
 \bfv{q}_{i} = \bfv{q}(i \DeltaT) \in \Omega, \; i \in \{ 0,\ldots,\LD \}
\end{equation}
and the fractional count matrix $\bfm{N}(\tau)$ of independent observed transitions by 
\begin{eqnarray}
	N_{ij}(\tau) & = & \frac{1}{\tauD}\sum\limits_{n = 0}^{\LD - \tauD} \chi_{i}\left(q_{n}\right) \chi_{j}\left(q_{n + \tauD}\right) 
\end{eqnarray}
%\begin{eqnarray}
%	C_{ij}(\tau, \beta) & = & \expect{\chi_i(t) \, \chi_j(t + \tau)}_{z(t, \beta)} \\
%	N_{ij}(\tau, \beta)	 & = & \int \limits_{0}^{L - \tau} \text{d}t \chi_i( z(t, \beta) ) \, \chi_j( z(t + \tau, \beta))
%\end{eqnarray}
%\begin{eqnarray}
%	N_{ij}(\tau, \beta)	 & = & \int \limits_{0}^{L - \tau} \text{d}t \chi_i( z(t, \beta) ) \, \chi_j( z(t + \tau, \beta)) \\
%	N_{ij}(\tau, \beta)	 & = & \sum\limits_{n = 1}^{L - 1} \chi_{i}(z_{n}(\beta)) \chi_{j}(z_{n + 1}(\beta)) 
%\end{eqnarray}
Assuming Markovianity, the likelihood that a given transition matrix produces the observations stored in the count matrix $\bfm{N}$ is
\begin{equation}
	\mathcal{L}(\bfm{T}) = \mathbb{P}(\bfm{N} | \bfm{T}) = \prod \limits_{i,j = 1}^{M} T_{ij}^{N_{ij}}.
\end{equation}
As a representative we choose the unique transition matrix $\hat{\bfm{T}}$ which maximizes this likelihood :
\begin{equation}
	\hat{T}_{ij} = \left[\text{argmax} \; \mathcal{L}(\bfm{T})\right]_{ij} = \frac{N_{ij}(\tau)}{\sum\limits_{k} N_{ik}(\tau)}.
\end{equation}

Alternatively, we can use the state-to-state time-correlation function $C_{ij}(\tau)$ given by
\begin{eqnarray}
	C_{ij}(\tau) &\equiv& \expect{\chi_i(0) \, \chi_j(\tau)}
\end{eqnarray}
which can be estimated in a similar fashion 
\begin{eqnarray}
	\hat{C}_{ij}(\tau) & = & \frac{1}{\LD - \tauD}\sum\limits_{n = 0}^{\LD - \tauD} \chi_{i}(q_{n}) \chi_{j}(q_{n + \tauD}) \\
	& = & \frac{\tauD}{\LD - \tauD} N_{ij}.
\end{eqnarray}
Although it can be formulated for different dynamical models\cite{Chodera:2009p6073} the present dynamical model 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 that sample from equilibrium, the correlation matrix will converge to a symmetric form $C_{ij}(\tau) = C_{ji}(\tau)$. In this case, we can use the estimator
\begin{equation}
	\hat{C}_{ij} = \frac{\tauD}{2 (\LD - \tauD)} (N_{ij} + N_{ji}) = \hat{C}_{ji}
	\label{eq:correlation-from-counts}
\end{equation}
and write the transition matrix in terms of the correlation matrix $\bfm{C}$ as
\begin{eqnarray}
 \hat{T}_{ij}(\tau) & = & \frac{\hat{C}_{ij}(\tau)}{\sum\limits_{k} \hat{C}_{ik}(\tau)}
 \label{eq:transitionmatrix-from-correlation}
\end{eqnarray}
which will also fulfill detailed balance.
%Given a partitioning of configuration space $\xi_{i}$ into $M$ subsets called states, lagtime $\tau$ and the temperature $\beta$ we need to find transition probabilities that maximize the probability to reproduce the observed state-to-state correlation functions, which can alternatively be regarded as a measure of uncorrelated transitions between states.

%Functions

%\begin{eqnarray}
%C_{ij}(\tau;\beta) & = & \frac{1}{Z} \int\limits_{z_0} \exp \left(  -\frac{H(z_0)}{k_B T} \right) \left( \chi_j(z_{\tau}) \chi_i(z_0)\right) \\
%C_{ij}(\tau) & \approx & \sum\limits_k^K \sum\limits_n^{N_k} w_{kn}(\beta) \cdot C_{ij}^{(kn)}(\tau) \\
%C_{ij}(\tau) & \approx & \sum\limits_k^K \sum\limits_n^{N_k} w_{kn}(\beta) \cdot \left\langle \chi_j(z_{kn}(t + \tau)) \chi_i(z_{kn}(t)) \right\rangle_{t} \\
%C_{ij}(\tau) & \approx & \sum\limits_k^K \sum\limits_n^{N_k} w_{kn}(\beta) \frac{1}{L-1} \sum\limits_{l=1}^{L-1} \chi_j(z_{kn,l+1}) \chi_i(z_{kn,l}) \\
%C_{ij}(\tau) & \approx & \sum\limits_k^K \sum\limits_n^{N_k} w_{kn}(\beta) \frac{1}{L-1} \sum\limits_{l=1}^{L-1} N_{ij}^{(kn)} \\
%C_{ij}(\tau) & \approx & \sum\limits_k^K \sum\limits_n^{N_k} w_{kn}(\beta) \hat{C}_{ij}^{(kn)}(\tau)
%\end{eqnarray}

\subsection{Transition probabilities from dynamical reweighting}

%\subsubsection{Dynamical reweighting (DR)}
\label{section:dynamical-reweighting}

We now demonstrate how transitions observed at all temperatures can be used to infer transition probabilities at any temperature of interest: For a canonical ensemble at inverse temperature $\beta = (k_{B} T)^{-1}$ the state-to-state correlation functions can be expressed as Boltzmann-weighted expectation functions
\begin{eqnarray}
C_{ij}(\tau;\beta) & = & \frac{1}{Z(\beta)} \int d \bfv{q}(0) d \bfv{p}(0) \exp \left(  -\beta \mathcal{H}(\bfv{q}(0),\bfv{p}(0)) \right) \left( \chi_i(\bfv{q}(0)) \chi_j(\bfv{q}(\tau)) \right)
\end{eqnarray}
where $Z(\beta)$ is the complete partition function of both kinetic and potential part.

Suppose we have a set of $N_k$ Hamiltonian trajectories $\z_{kn}(t)$, $n = 1,\ldots,N_k$, $t \in [0,T]$, in which the initial phase space points $\z_{kn}(0)$ are sampled from canonical (NVT) distributions at corresponding inverse temperatures $\beta_k$, $k = 1,\ldots,K$. By the application of \emph{dynamical reweighting}\cite{Shirts:2008p2350}~, a correlation function $C_{ij}(\tau;\beta)$ can be estimated using the entire set of trajectories at all temperatures as
\begin{eqnarray}
	\hat{C}_{ij}(\tau; \beta) & \approx & \sum\limits_k^K \sum\limits_n^{N_k} w_{kn}(\beta) \cdot \hat{C}_{ij}^{(kn)}(\tau)
\end{eqnarray}
where the individual microcanonical contributions to the correlation functions are given by Eq. (\ref{eq:correlation-from-counts})
\begin{equation}
	\hat{C}_{ij}^{(kn)} = \frac{\tauD}{2 (\LD - \tauD )} \left( N_{ij}^{kn} + N_{ji}^{nk} \right)
	\label{eq:reversibleCounting}
\end{equation}

%\begin{eqnarray}
%	C_{AB}(\tau;\beta) &\equiv& [Z(\beta)]^{-1} \int d\z_0 \, e^{-\beta H(\z_0)} \, \mathcal{A}[\z(t)] \nonumber \\
%&\approx& \sum\limits_{k=1}^K \sum\limits_{n=1}^{N_k} w_{kn}(\beta) \, \mathcal{A}[\z_{kn}(t)]
%\end{eqnarray}
%where $\mathcal{A}[\z(t)]$ is a functional of the trajectory $\z(t)$ whose expectation produces the desired correlation function, such as $\mathcal{A}[\z(t)] = A(\z(0)) \, B^*(\z(\tau))$.

%To estimate the transition probabilities $T_{ij}(\tau;\beta)$, we first compute an estimate of the correlation functions $C_{ij}(\tau;\beta)$
%\begin{eqnarray}
%\hat{C}_{ij}(\tau;\beta) &=& \sum_{k=1}^K \sum_{n=1}^{N_k} w_{kn}(\beta) \, A_{kn}(\tau)
%\end{eqnarray}
%where we have defined the per-trajectory observable $A_{kn}(\tau)$ to make use of overlapping trajectory segments (in the case $T > \tau$) and the time-reversibility of Hamiltonian dynamics, thus ensuring $\hat{C}_{ij}(\tau;\beta) = \hat{C}_{ji}(\tau;\beta)$:
%\begin{eqnarray}
%A_{kn}(t) &\equiv& \frac{1}{2(T-t)} \int_0^{\tau-t} dt_0 \, [ \chi_i(\z_{kn}(t_0)) \, \chi_j(\z_{kn}(t_0+t)) \nonumber \\
%&&\mbox{} + \chi_j(\z_{kn}(t_0)) \chi_i(\z_{kn}(t_0+t)) ]
%\end{eqnarray}
The row-stochastic transition matrix estimate $\hat{\bfm{T}}(\tau;\beta)$ is then estimated from (\ref{eq:transitionmatrix-from-correlation})
\begin{eqnarray}
\hat{T}_{ij}(\tau;\beta) &=& \frac{\hat{C}_{ij}(\tau;\beta)}{\sum\limits_{k} \hat{C}_{ik}(\tau;\beta)}
\end{eqnarray}
As $\hat{\bfm{C}}(\tau;\beta)$ is symmetric by construction, $\hat{\bfm{T}}(\tau;\beta)$ will be reversible (satisfy detailed balance).

The normalized trajectory weights $w_{kn}(\beta)$ are specified in terms of unnormalized weights, which can be computed by
\begin{eqnarray}
	w_{kn}(\beta) & = & \tilde{w}_{kn}(\beta) / \sum\limits_{k=1}^K \sum\limits_{n=1}^{N_k} \tilde{w}_{kn}(\beta) \\
	\label{eq:normalized-trajectory-weights}
	\tilde{w}_{kn}(\beta) &=& \left[ \sum_{k'=1}^K N_{k'} \, \exp[\hat{f}_{k'} - (\beta - \beta_{k'}) \, E_{kn}] \right]^{-1}
	\label{eq:unnormalized-trajectory-weights}
\end{eqnarray}
where $E_{kn} \equiv H(\z_{kn}(0))$ denotes the total energy of the trajectory, which is constant over trajectories for Hamiltonian dynamics\cite{Chodera:2009p6073}. The dimensionless free energies $\hat{f}_i = - \ln Z_i + c$ are determined by the solution of a set of self-consistent equations
\begin{eqnarray}
	\hat{f}_i &=& - \ln \sum_{k=1}^K \sum_{n=1}^{N_k} \tilde{w}_{kn}(\beta_{i}), \; \forall \; i \in \{ 1, \ldots, K \}
	\label{eq:free-energies}
\end{eqnarray}
which can be obtained efficiently by a number of means~\cite{Shirts:2008p2350}, although it is often necessary to work with logarithmic representations to avoid numerical instability. A complete derivation of all expressions above is presented in Ref. \cite{Chodera:2009p6073}.

%{\color{red}[JHP:] CITE THEORY PAPER OR COPY APPENDIX	}

\subsection{Estimation of uncertainties in transition probabilities}
\label{appendix:propagation-of-error}

The statistical uncertainty in $\hat{C}_{ab}(\tau,\beta)$ can be estimated in a straightforward manner using the machinery of MBAR~\cite{Shirts:2008p2350}. We following the derivation presented in \cite{Chodera:2009p6073} and first reindex $k$ and $n$ in the datasets $\{\z_{kn}, E_{kn}, \hat{C}_{ab,kn}\}$ by replacing it with one single index $n = 1,\ldots,N$, where $N = \sum_{k=1}^K N_k$ and $a,b = 1,\ldots,M$ with $M=6$ being the number of states in the model considered. This results in the datasets $\{\z_n, E_n, \hat{C}_{ab,n}\}$.

We form the $N \times (K+1+M^2)$ matrix $W_{nk}$ such that
\begin{eqnarray}
W_{nk} = \begin{cases}
\exp[-\beta_k E_n] \quad & k = 1,\ldots,K \\
\exp[-\beta E_n] \quad & k = K+1 \\
\hat{C}_{ab,n} \, \exp[-\beta E_n] \quad & \parbox{2.0in}{$k = K+1+a(M-1)+b$, \\ $a,b = 1,\ldots,M$}
\end{cases}
\end{eqnarray}
The uncertainty of the transition probabilities $T_{ij}(\hat{C}_{11}, \ldots, \hat{C}_{MM})$ can then be estimated by
\begin{eqnarray}
\delta^2 T_{ij} &\equiv& \expect{(T_{ij} - \expect{T_{ij}}_\beta)^2}_\beta \nonumber \\
&\approx& \sum_{a,a^{'},b,b^{'}=1}^M \left[ \frac{\partial T_{ij}}{\partial \hat{C}_{ab}} \right] \, \left[ \frac{\partial T_{ij}}{\partial \hat{C}_{a'b'}} \right] \, \delta \hat{C}_{ab} \delta \hat{C}_{a'b'} 
\end{eqnarray}
where the covariance of the estimates $\hat{C}_{ab}$ and $\hat{C}_{a'b'}$ can be estimated as
\begin{eqnarray}
\delta \hat{C}_{ab} \delta \hat{C}_{a'b'} &\approx& \hat{C}_{ab} \, \hat{C}_{a'b'} [ \hat{\Theta}_{K+1,K+1}  + \hat{\Theta}_{K+1+a(M-1)+b,K+1+a'(M-1)+b'} \nonumber \\	
&&\mbox{} - \hat{\Theta}_{K+1,K+1+a'(M-1)+b'} - \hat{\Theta}_{K+1+a(M-1)+b,K+1}] \label{equation:covariance-of-estimates}
\end{eqnarray}
with covariance matrix $\hat{\bfm{\Theta}} = [ (\bfm{W}^\T \bfm{W})^+ - \bfm{N} ]^{+}$ and $\bfm{N} = \mathrm{diag}(N_1, \ldots, N_K, 0, \ldots, 0)$. The~${}^{+}$ here denotes the generalized inverse. Using Eq. (\ref{equation:transition-matrix-definition}) the sensitivity of $T_{ij}$ to the correlation $\hat{C}_{ab}$ is given by
\begin{equation}
\frac{\partial 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},
\end{equation}
with $\hat{C_{i}} = \sum_{i} C_{ij}$ being the row sum of the correlation matrix $\bfm{C}$.
The final expression for the covariance matrix of transition probabilities $\delta^2 T_{ij}$ is complex but still enumerable. A more thorough derivation of the above can be found in \cite{Chodera:2009p6073,Shirts:2008p2350}.

\subsection{Modified parallel tempering protocol}
\label{section:modified-parallel-tempering-protocol}

By using a modified parallel tempering protocol, a set can be generated of Hamiltonian trajectory segments $\z_{kn}(t)$ of uniform length $T \ge \tau$ the initial phase space points $\z_{kn}(0)$ of which are sampled from the canonical (NVT) ensemble at corresponding inverse temperatures $\beta_1,\ldots,\beta_K$.

We start by assuming that some process was used to generate the initial phase space points $\z_{k0}(0)$ from equilibrium at each corresponding inverse temperature $\beta_k$
\begin{eqnarray}
\z_{k0}(0) &\sim& [Z(\beta_k)]^{-1} \, e^{-\beta \mathcal{H}(\z_{k0}(0))}
\end{eqnarray}
This initial phase space point generation may be by means, for example, of a standard parallel tempering protocol, or by running the modified protocol for a number of iterations starting from a single configuration.

Consider iteration $n$ of the algorithm. For each temperature index $k = 1,\ldots,K$, Hamilton's equations of motion are propagated using a symplectic integrator to generate trajectories of $\z_{kn}(t)$ of length $T$. Finally, we propose exchanges\cite{Chodera:2006p2352} between the final configurations $\z_{in}(T)$ and $\z_{jn}(T)$ 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. The Metropolis-like probability\cite{Sugita_CPL314_141} of accepting or rejecting the exchange depends on the final potential energies of the configurations $U_i$ and $U_j$ with
\begin{eqnarray}
	P_\mathrm{exch}(U_i, \beta_i; U_j, \beta_j) &=& \min \left\{ 1, \exp[-(\beta_i - \beta_j) (U_j - U_i)] \right\} \nonumber
\end{eqnarray}
Regardless of whether the exchange is accepted or rejected, we reassign the velocities according to the Maxwell-Boltzmann distribution at the new (old) 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:modified-parallel-tempering-protocol-proof}).
This satisfies the conditions defined by Okamoto\cite{Sugita_CPL314_141} in order for the kinetic energies to not appear in $P_\mathrm{exch}$. The reason for reassignment of velocities instead of rescaling is that, when using Hamiltonian trajectories no thermostating would otherwise take place.

\subsection{Bayesian estimate of transition probabilities from a single temperature}

To analyze how the combination of data from all temperatures reduces the uncertainty in the estimation of the transition matrix estimates we compare the present method to two Bayesian methods of error estimation working on data only at a single temperature. 
Both Bayesian methods sample transition probabilities or rates according to the same likelihood function, but employ 
different proposal steps and, more importantly different prior probability distributions.

%\subsubsection{Transition Matrix Estimation (TE)}
\label{section:transition-matrix-estimation}
Let $N_{ij}(\tau,\beta)$ be the transition counts from all trajectories at inverse temperature $\beta$. Then the a posteriory probability of transition matrix $T_{ij}$ given this observation is
\begin{equation}
\mathbb{P}(\bfm{T}|\bfm{N}) \propto \mathbb{P}(\bfm{T}) \; \mathbb{P}(\bfm{N}|\bfm{T})  = \mathbb{P}(\bfm{T}) \prod_{i,j \in S} T_{ij}^{N_{ij}}.
\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}(T) \equiv \prod_{i,j \in S} 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}.
The approach proposed in Ref. \cite{Buchete:2008p2354} does not estimate the transition probabilities directly, but uses a likelihood function based on parameters of a reversible rate matrix $\bfm{K}$ with the corresponding likelihood
\begin{equation}
	\mathbb{P}(\bfm{K}|\bfm{N}) \propto \mathbb{P}(\bfm{K}) \prod_{i,j \in S} \exp \left( \tau \bfm{K}\right) _{ij}^{N_{ij}}
	\label{eq:hummer-likelihood}
\end{equation}
and the detailed balance constraint
\begin{equation}
	\pi_{i} K_{ij} = \pi_{j} K_{ji}.
\end{equation}
The choice of parameterization for $\bfm{K}$ in Ref. \cite{Buchete:2008p2354} stores the logarithms of the upper-right triangular matrix (without diagonal entries) and the equilibrium distribution, thus assuring a rate matrix with negative eigenvalues and also positive off-diagonal rates. After sampling rate matrices with a Metropolis-Monte Carlo scheme the related set of transition matrices with the lagtime $\tau$ is computed  by
\begin{equation}
	\bfm{T} = \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} 
%{\color{red} Do we really need this?}.
\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 302K \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}
To illustrate the dynamical reweighting method described above, we apply it to estimate transition probabilities between conformational states for the terminally-blocked alanine peptide (Ace-Ala-Nme). 

For the analysis a parallel tempering dataset for the terminally-blocked alanine peptide in a box of 431 TIP-3P water molecules was generated as described previously \cite{Chodera:2006p2352}. The parallel tempering simulation was conducted according to the modified protocol described in the theory section \ref{section:modified-parallel-tempering-protocol} generating an ensemble of 501 Hamiltonian trajectories of 20 ps length at each of 40 temperatures spanning 273 to 600 K, with peptide configurations stored every 0.1 ps. The temperatures were exponentially spaced, thus providing good overlap in the potential and total energy distributions between neighboring temperatures and resulting in an average exchange acceptance probability of $\sim50\%$ (see Figure~\ref{figure:energydistribution}). 
\begin{widetext}
\begin{figure}[]
\includegraphics[width=16cm]{fig/plot_EnergyOverlap_InitialEnergy_Gauss.pdf}
\caption{\label{figure:energydistribution} Distribution of energies in the trajectories in $kJ/\text{mol}$ for all 40 temperatures. Highlighted are temperatures 302~K, 426~K and 600~K to allow reference with convergence properties in Figure \ref{fig:convergence-of-bayesian-methods}. The temperature spacing was chosen as exponential so as to provide an average overlap of $\sim50\%$ between neighboring temperatures.}
\end{figure}
\end{widetext} 
A velocity Verlet integrator\cite{swope:1982a} (with bonds involving hydrogen atoms constrained) was used to produce the dynamical trajectories. The fluctuations in total energy over the 20 ps trajectories was minimal and the drift negligible (see Table \ref{table:pt-energy-rms-and-drift}). The production run followed a 1 ns equilibration phase, ensuring that all initial configurations were drawn from equilibrium at their respective temperatures.
\begin{table}
\begin{ruledtabular}
\begin{tabular}{lrrrr}
Temp [K] & RMS [J/mol] & rel. RMS & Drift [J/(ps mol)] & rel. Drift [1/ps]\\ \hline
302 & $214.844 \pm 0.648$ & $6.41 \cdot 10^{-5}$ & $0.608  \pm 0.249 $ & $1.82 \cdot 10^{-7}$  \\
600 & $26.617 \pm 0.083$ & $1.67 \cdot 10^{-5}$ & $1.061 \pm 0.432 $ & $6.66 \cdot 10^{-7}$
\end{tabular}
\end{ruledtabular}
\caption{\label{table:pt-energy-rms-and-drift}RMS and Drift of Energies over 20ps trajectories}
\end{table}
Previous work has demonstrated that a Markov model based on a six-state decomposition as depicted in Figure~\ref{fig:energyLandscape} can accurately describe the dynamics of this peptide for lagtimes longer than $\tau=6 \text{ ps}$\cite{Chodera:2006p2352}. In order to conveniently directly apply the reweighting procedure we used the same state decomposition for all temperatures. 

To evaluate the accuracy of the methods for estimating transition probabilities, we compare the separate estimates from (DR),(TE) and (RE) to a reference simulation of 6x10000 short, 10~ps trajectories shot out of the equilibrium of each state at 302 K. The PT simulation in comparison furnishes a total of 1500 independent transitions shared between all states at that temperature, which means that the shooting dataset has 40 times more observations compared to the PT simulation, thus providing a statistically reliable reference.

The system is small enough such that a reasonable statistics can be obtained with moderate CPU requirements while at the same time being complex enough such that some transitions (and even some states) are sampled only at high temperatures. 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 a single temperature only.
\begin{figure}[]
\includegraphics[width=8cm]{fig/plot_energyLandscapeAt302K.pdf}
\caption{\label{fig:energyLandscape} Terminally-blocked alanine potential of mean force and state boundaries. 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 parallel tempering simulation using WHAM\cite{Chodera:2007p2346,Kumar:1992p860}, truncated at $10 k_{B}T$ (white regions). The six manually identified states are labelled in black\cite{Chodera:2006p2352}.}
\end{figure}

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

A comparison of the transition probabilities between all 6x6 pairs of states as a function of temperature is given in Figure~\ref{fig:transition-probabilities-temp}. The blue solid lines give the estimates from dynamical reweighting (DR)\cite{Chodera:2006p2352} using all available data at all temperatures as described in the Theory section.
To obtain the dimensionless free energy estimates $\hat{f}_{i}$ we solved the set of self-consistent equations in (\ref{eq:free-energies}) with a relative convergence tolerance in the residual of $10^{-7}$. Transition probabilities were also estimated at one intermediate temperature between each pair of simulated temperatures. 

The red dashed lines in Figure~\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), depicted by dotted green lines, the sampling as proposed in Ref. \cite{Buchete:2008p2354} was used. After an equilibration phase again a total of 10,000 samples were stored for each temperature separately. For the convergence of both methods see supplementary Figure~\ref{fig:convergence-of-bayesian-methods}. The black cross-hair in Figure~\ref{fig:transition-probabilities-temp} refers to the reference values at 302 K. 
Qualitatively all methods agree, especially in the lower energy states (1 to 4). However, the reweighting estimate, which uses the combined data from all temperatures, has smaller uncertainties than the estimators that use only individual temperatures.
The overall agreement with the reference simulation is best for dynamical reweighting. The two Bayesian methods have almost identical predictions for states with many observed transitions, which is as expected since they use the same likelihood functions and the influence of the choice of the prior probability distribution is minimal here. However, the differences in the states with few observations (states 5 and 6) arise from the influence of different prior probability distributions.

\begin{widetext}
\begin{figure}[!htb]
\includegraphics[width=\textwidth]{fig/plot_allTransitionProbabilites_woLimits.pdf}
\caption{\label{fig:transition-probabilities-temp} Comparison of all 6x6 inter-state transition probabilities as a function of temperature with error bars showing $95\%$ confidence intervals. Blue lines show the transition probabilities estimated using the present method (DR). Red lines show the estimates from transition matrix estimation (TE) 
%, the green lines resp. the rate matrix estimation (RE) both 
computed from only single temperature data. 
The black cross-hair indicates the reference using the shooting trajectory data (ST) at 302K only. It is apparent that the reweighting method provides useful and bounded estimates across the whole temperature range even at temperatures where the corresponding transition was not observed (e.g., at transitions from state 5 and 6 for low temperatures). The Bayesian counting method without reweighting is very noisy and not useful in situations at low temperatures for transitions with poor statistics.}
\end{figure}
\end{widetext}

%Surprisingly, there is still a lot of information about transitions to/from states that haven't been sampled (or have been poorly sampled) from the Bayesian Sampling approach at single temperatures, due to the incorporation of the reversibility constraint. 

%For transitions, that are not sampled at certain temperature ranges, the maximum-likelihood estimate obtained with our reweighting method is approximately zero (s) 
%In these cases the uncertainty is underestimated, since it is derived from the low temperatures for poorly sampled transitions (e.g. 3->6), suspect reweighting may significantly underestimate uncertainty since it is derived from the local (hessian) and thus for values close to zero a bad approximation for the standard deviation of this highly unsymmetric distribution (s. fig ...)

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

For a detailed comparison the Bayesian analysis method with reversibility constraint for transition matrices\cite{Noe:2008p2459} (TE) was applied to the set of shooting trajectories at 302 K as well to generate an estimate of reference transition probabilities obeying detailed balance. The results of the comparison at 302 K between the different estimation methods are shown in Figure~\ref{fig:detailed-probabilities-302K}. All colors are the same as in Figure~\ref{fig:transition-probabilities-temp}.

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 Figure~\ref{fig:transition-probabilities-temp}). Generally, for transition probabilities close to zero or one the normal distribution is a poor approximation to this highly asymmetric distribution and therefore leads to too large uncertainties (see Figure~\ref{fig:detailed-probabilities-302K}) in cases. 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.

The reweighting method performs overall very well compared to the single-temperature estimates. Even transition probabilities that are sampled very poorly at 302~K (such as for transitions involving states 5 and 6) have a 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). The dynamical reweighting has a smaller deviation than both Bayesian methods for both high- and low-energy states.
\begin{table}
\begin{ruledtabular}
\begin{tabular}{lccc}
 & \multicolumn{3}{c}{Variance in Estimation} \\ \cline{2-4}
 & (DR) & (TE) & (RE)\\ \hline
Low-energy states (1 to 4) & 0.007 & 0.023 & 0.055 \\
High-energy states (5 and 6) & 0.140 & 0.267 & 0.288 \\ \hline
All Transitions & 0.079 & 0.151 & 0.167
\end{tabular}
\end{ruledtabular}
\caption{\label{table:comparison-302K} Standard Deviation in the absolute difference of transition probabilities compared to the reference simulation (ST) at 302K for the three methods of Markov model estimation and high- and low-energy subsets of transitions. Dynamical reweighting (DR) shows the smallest spread in deviation from the reference simulation (ST).}
\end{table}

\begin{figure}[]
\includegraphics[width=\textwidth]{fig/plot_allComparisonDetailed302K_woLimits.pdf}
\caption{\label{fig:detailed-probabilities-302K} Detailed Comparison of transition probabilities at 302 K. 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}

%Any observations about metastable state structure or grouping as temperature changes? (PCCA+ as a function of temperature?) [JHP: We could try that, but I don?t think there is a way of visualizing this. Maybe show the shift of barriers between metastable states with temperature? We can only do this with the 6 state model, since all dynamical reweighting is based on this. We could however plot the stability vs temperature and show, that these decrease with temperature, but this is already included in the diagonal elements of the transition matrix plot! So in a way redundant]

%The blue lines show the estimated transition probabilities using the present method, the result. Error bars indicate the 95\% confidence intervals. For comparison the black cross-hair shows the maximum-likelihood transition probabilities computed from the shooting trajectories. It is apparent that the reweighting method (DR) provides useful and bounded estimated across the whole temperature range even when the corresponding transition was not observed at that temperature. Both Bayesian methods without reweighting are very noisy and provide worse estimates for all transitions especially at lower temperatures.

%\begin{figure}[]
%\includegraphics[width=8cm]{fig/plot_DeviationInTransitionProbability_All.pdf}
%\caption{\label{fig:comparison-302K} Relative Difference of Transition Probabilities to the reference values at 302K for each type of Markov model estimation. Blue: dynamical reweighting (DR), Red: transition matrix estimation (TE), Green: rate matrix estimation (RE). The dynamical reweighting show the smallest spread in deviation from the reference simulation (ST).}
%\end{figure}

\subsection{Comparison of temperature dependence of eigenvalues}
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}^{*}$ indicated by the corresponding eigenvector by
\begin{equation}
	t_{i}^{*} = - \tau / \log(\lambda_{i}).
\end{equation}
Hence, eigenvalues close to one imply slow processes i.e., those we are mostly interested in. We investigated the dependence of the eigenvalues on the temperature in the present system. Figure~\ref{fig:eigenvalues-temp} compares estimates for the second and third eigenvalue $\lambda_{2}, \lambda_{3}$ of the transition matrix estimated 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. To estimate the errors of DR we used linear error propagation of the uncertainties in the transition matrix to the errors in the eigenvalues \cite{Hinrichs:2007p5942}.
\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} Red: Single temperature estimation of transition matrix (TE), Blue: dynamical reweighting estimation (DR), Black: reference using shooting trajectories (ST). Left: Comparison of the second eigenvalue vs temperature, Right: Comparison for third eigenvalue. The third eigenvalue is well predicted by both estimation methods at all temperatures, while the second eigenvalue at low temperatures is only detected by dynamical reweighting.}
\end{figure}
At low temperatures, 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 the good agreement of the expectation values with the reference. We assume that the inappropriate approximation of the asymmetric distributions with normal distributions used for the linear error propagation lead here as well to an overestimation of the errors in the transition probabilities.
%{\color{red} [JHP:] Not completely sure, where these really come from? But coincide very well with "gold standard" at 302 K. Is this coincidence, since the error bars are so large? Maybe we just mention this and not give an explanation. I played a little with exact eigenvalue distributions from dirichlet distributed transition matrices. This show, that the error in the gauss approximation is especially overestimated for values close to one.}

The third largest eigenvalue is predicted by both methods equally well (Fig.~\ref{fig:eigenvalues-temp}). Although this was found as the second-largest eigenvalues in the single-temperature estimates, which missed sampling the slowest process 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 mixed up for faster processes. The comparison of Markov Models is a nontrivial task, for which we used a symmetrized form of the transition matrix $\bfm{T}^{sym}$ and expand it into a sum of rank one matrices $\bfm{Q}_{i}$ 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} \\
   	              & = & \sum \limits_{i=1}^{M} \bfm{Q}_{i}.
\end{eqnarray}
Here, $\bfv{\pi} = \{ \pi_{1}, \ldots, \pi_{M}\}$ is the equilibrium distribution, $\bfm{R} = \{ \bfv{r_{1}}, \ldots, \bfv{r_{M}} \}$ the matrix of eigenvectors of $\bfm{T}^{sym}$ and $\lambda_{i}$ the corresponding eigenvalues, which are equal to the eigenvalues of $\bfm{T}$. Each of the submatrices $\bfm{Q}_{i}$ can be considered as a part of the full transition matrix working on a timescale given by the respective eigenvalue. These subprocesses are similar if their eigenvectors match, as can be determined by the scalar product. 
\begin{widetext}
\begin{figure}[]
\includegraphics[width=\textwidth]{fig/plot_Eigenvector_Matching_Comparison.pdf}
\caption{\label{fig:eigenvectorMatchingComparison} Similarity matrices (scalar product) of eigenvectors from symmetrized transition matrices estimated with different methods (ST, DR, TE, RE) at 302 K. The eigenvectors indicate the states involved in the process, thus high similarity (red) indicates a good approximation to the reference process (ST). Eigenvectors are sorted as descending eigenvalues. The 2nd eigenvector is found correctly only by dynamical reweighting, meaning that the single temperature estimations are unable to correctly predict the slowest process at 302~K.}
\end{figure}
\end{widetext}

\subsection{Contributions from different temperatures to the estimates of expectation values}

The contribution of each trajectory snippet to the estimation of any expectation value at any given temperature is illustrated in Fig.~\ref{fig:contribution-to-weights}. The left plot shows the average normalized weights 
\begin{equation}
	\bar{w}_{kl} = \frac{1}{N_{k}} \sum\limits_{n} w_{kn}(\beta_{l})
	\label{eg:average_weight}
\end{equation} 
given in (\ref{eq:unnormalized-trajectory-weights}), where, on average, seven temperatures contribute more than 1\% to the expectation. The right hand plot illustrates the contribution to the transition counts $N_{65}$ for the transition $6 \rightarrow 5$ given by
\begin{equation}
	\frac{1}{N_{k}} \sum \limits_{n=1}^{N_{k}} \tilde{w}_{kn}(\beta_{l}) \hat{C}_{ij}^{kn},
\end{equation}
with $k$ indicating the temperature contributed from, $l$ the temperature estimated at and $i,j$ representing the transition $i \rightarrow j$.
As expected most information contributing to a specific temperature estimate is always contained in the simulations performed at the nearest temperatures. 
\begin{widetext}
\begin{figure}[]
\includegraphics[width=5cm]{fig/plot_normalized_weights_Wkn.pdf}
%\includegraphics[width=5cm]{fig/plot_contribution_to_summary.pdf}
\includegraphics[width=5cm]{fig/plot_contribution_to_estimate_W65kl.pdf}
\caption{\label{fig:contribution-to-weights} Left: 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}$  averaged over all trajectories at the same temperature. On average, 7 temperatures contribute more than 1\% each to the estimation. Right: Contributions to the estimation of transition counts for the transition $6\rightarrow 5$. The sum of one row equals to the total counts estimated by the method at the desired temperature.}
\end{figure}
\end{widetext}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DISCUSSION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Discussion}

%\subsection{Assessment of reweighting}
The present method provides a means of generating an estimate of transition probabilities from parallel tempering MD simulations of biomolecules as a continuous function of temperature. Even at intermediate temperatures not simulated at, the error bars are much smaller then obtained with either single-temperature methods. At low temperatures, where some transitions are not observed at all, non-zero transition probabilities can still be estimated. %although the method tends to underestimate probability of transitions with very low statistics.

Additionally, the method allows to differentiation of the estimates of transition probabilities with respect to the inverse temperature $\beta$, because the trajectory weights $w_{kn}(\beta)$ are differentiable functions with respect to  temperature (\ref{eq:unnormalized-trajectory-weights}). This allows, in principle, thermodynamic properties to be computed (e.g., heat capacities), although these quantities are in general numerically difficult to handle, since the trajectory weights can easily span hundreds of orders of magnitude. 

The method requires a set of parallel tempering simulations with a modified protocol to produce a series of NVE trajectories with initial configurations drawn from the NVT ensemble. For very large systems, the PT simulation might not be long enough to globally converge. In this case we cannot use reversible counting as in Eq. (\ref{eq:reversibleCounting}) to enhanced the statistics, but the method can still be applied without the detailed balance constraint as long as we draw from equilibrium inside each set $\Gamma_{i}$.
The PT requirement of good exchange rates also ensures good overlap in the contribution to the dynamical reweighting for neighboring temperatures.  
The approach itself is not limited to Hamiltonian trajectories but can be extended to other dynamics as long as an analytical connection between the weighting factors and the temperature exists (e.g. Boltzmann distributed probabilities)\cite{Chodera:2009p6073}. For example, this is the case for Brownian and Langevin dynamics. 

%\subsection{Thermally-activated processes and entropic barriers}

The degree to which the use of parallel tempering can speed up convergence is a matter of discussion. 
%{\color{red} [JCS] Needs: References} 
While activated processes will be sampled more often at higher temperatures, entropic barriers become less and less probable to pass, effectively limiting the possible increase in simulation speed. The same problem limits the range of contributing temperatures here, too. Nevertheless, information about activated processes is transferred from higher to lower temperatures and, for entropic barriers from lower to higher temperatures. Transitions in the alanine peptide are dominated by activated processes (i.e., enthalpic barriers).

%\subsection{Comments on performance of Bayesian estimate}

Both single-temperature methods give quite similar results for transitions with good statistics, differing mostly for transitions that have only rarely been sampled due presumably to the influence of different prior probability distributions. Surprisingly, the Bayesian estimates provide a reasonable bound on transition probabilities to and from a state even when the state is not even sampled. 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.

There is, however, a dependence on the prior, which leads to different predictions in cases where the state is not or only rarely sampled. The rate matrix estimation (RE) assumes, in addition to the detailed balance constraint, positivity of all eigenvalues and non-negative off-diagonal entries. The uniform distribution of parameters in logarithmic space leads most likely to favoring of small transition probabilities in states with poor transition statistics.

%\subsection{Requirements for convergence}

The way in which the transition probabilities are estimated in terms of equilibrium correlation functions requires that the trajectory segments sampled during the parallel tempering simulation are drawn from the equilibrium distribution and that the trajectories to be reweighted are uncorrelated so as to permit an estimation of the statistical error.

%While that is the case here (see Figure~\ref{fig:convergence-of-bayesian-methods}), this requirement may be more stringent than is necessary for transition probabilities.

%{\color{red} Bayesian method requires only ``local'' equilibrium within each state; we could presumably modify the reweighting scheme to only require this as well. [JDC: For example, we could select all of the trajectories that originate or terminate within a given state i and, from these, apply dynamical reweighting to estimate $\left\langle \chi_{j}(x_{t}) \right\rangle$ over this set of trajectories. The problem is that we would lose the built-in detailed balance constraint, unless we also enforced $p_{i} T_{ij} = p_{j} T_{ji}$ by ``the method of control variates'', which increases the complexity of the estimation procedure significantly. It may be best to just mention this idea and leave it alone, and go on to the Bayesian reweighting procedure instead.]}

%\subsection{Outlook}

The predictions of mean values are very good, while the quality of the error estimation is limited to a Gaussian approximation, which is problematic in cases where transition probabilities are close to 0 or 1 and the probability distribution is thus very asymmetric. 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.
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.
JHP and FN acknowledge support from DFG Research Center Matheon. JCS acknowledges funding from the U.S. Department of Energy ``Multiscale Mathematics'' {\color{red} Grant No. missing!}

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

\bibliographystyle{prsty} \bibliography{markov-model-from-parallel-tempering,markov-model-from-parallel-tempering-prinz-all}

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

\appendix
\begin{widetext}

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

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

Here we prove that the modified parallel tempering protocol described in the theory section samples from the canonical stationary distribution at all temperatures.

We define stationary distributions for momenta $\p$ and coordinates $\q$ in Cartesian space $\Re^{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 the current phase space points of which 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 \\
&+& \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)} - \min\left\{ \frac{e^{-\beta_1 U(\q'_1)}}{Q(\beta_1)} , \frac{e^{-\beta_2 U(\q'_1)}}{Q(\beta_1)} \frac{Q(\beta_1)}{Q(\beta_2)} \right\} + \min\left\{ \frac{e^{-\beta_2 U(\q'_1)}}{Q(\beta_2)} , \frac{e^{-\beta_1 U(\q'_1)}}{Q(\beta_2)} \frac{Q(\beta_2)}{Q(\beta_1)} \right\} \nonumber \\
&=& \pi_q(\q'_1 | \beta_1)  - \min\{ \pi_q(\q'_1 | \beta_1), \pi_q(\q'_1 | \beta_2)\} + \min\{ \pi_q(\q'_1 | \beta_2, \pi_q(\q'_1 | \beta_1) \nonumber \\
&=& \pi_q(\q'_1 | \beta_1)
\end{eqnarray}
Therefore, after the exchange attempt, the new configuration $\q'_1$ is still at equilibrium at the inverse temperature $\beta_1$.
(A similar series of steps can be applied for the temperature $\beta_2$.)

Again redrawing the momentum from a 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 figure \ref{fig:convergence-of-bayesian-methods} for various temperatures. Both methods converge after approximately 5000 samples.

\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: 273K, Yellow: 426K, Red: 600K. After about 5000 samples the confidence intervals stabilize suggesting reasonably well sampled transition probabilities.}
\end{figure}

\clearpage

%\section{Things to clarify (NOT IN THE PAPER !!!)}
%[Explain: What is meant with additive effects of uncorrelated statistical error?] [JDC: If I draw independent samples from two distributions $\pi_{1}(x)$ and $\pi_{2}(x)$ and try to estimate $E_{1}[A] - E_{2}[A]$ by estimating $E_{1}[A]$ using the samples from $\pi_{1}$ only and $E_{2}[A]$ using the samples from $\pi_{2}$ only, the squared-uncertainty in my estimates of $A$ for $\pi_{1}$ and $\pi_{2}$ will be additive. If I instead use all the data to estimate $E_{1}[A] - E_{2}[A]$ using reweighting, this error can be reduced, in some cases drastically.]



%\section{Simulation Setup (from MMS)}

%Using the {\tt LEaP} program from the AMBER7\cite{AMBER7} molecular mechanics package, a terminally blocked alanine peptide(sequence ACE-ALA-NME; see Figure 1)was generated in the extended conformation with peptide force field parameters taken from the AMBERparm96 parameter set\cite{AMBER-parm96}. The system was subsequently solvated with 431TIP3P water molecules\cite{jorgensen:1983a} in a truncated  octahedral simulation box with dimensions chosen to ensure all box boundaries were at least $7\text{\AA}$ from any atom of the extended peptide. All minimization and molecular dynamics simulations were conducted using the sander program from the  AMBER7 package. Default nonbonded cutoffs were used, bonds to hydrogen were constrained with SHAKE using a tolerance of $10^{-8}$ \cite{SHAKE}, and long-range electrostatics were treated by the particle-mesh Ewald method \cite{ewald93} with the default settings. 

%The system was first subjected to 50 steps of steepest descent energy minimization, followed by 1000 steps of conjugate gradient optimization. To equilibrate the explicit solvent system to the appropriate volume, a 100 ps molecular dynamics simulation was performed with the temperature adjusted to 300 K and the pressure to 1 atm by the Berendsen weak-coupling algorithm \cite{berendsen:1984a} with temperature and pressure relaxation time constants of 1 ps and 0.2 ps, respectively. The simulation box was fixed at the final size obtained from this equilibration step, with a volume of 13 232 $\text{\AA}^3$, in all subsequent simulations.

%\section{Parallel Tempering Protocol (from MMS)}

%In order to broadly explore the configuration space of the peptide and ensure that all important conformational substates were located, a parallel tempering (or replica exchange among temperatures) molecular dynamics simulation \cite{Sugita_CPL314_141} was conducted using a parallel Perl wrapper for the sander program \cite{code:rex}. Forty replicas were used, with replica temperatures exponentially distributed over the range $273-600 \text{K}$, yielding an average exchange acceptance probability of about $50\%$. All momenta were reassigned from the Maxwell-Boltzmann distribution at the appropriate replica temperature after each exchange attempt, and constant-energy, constant-volume molecular dynamics with a 2 fs timestep was performed between exchange attempts. The algorithm used to select pairs of replicas for temperature exchange attempts starts from the highest temperature replica, attempts to swap the configuration for the next-lowest temperature replica using the Metropolis-like criteria, and proceeds down the temperatures in this manner. On the next iteration, swapping attempts start from the lowest temperature and proceed 
%upward, and this alternation in direction is continued in subsequent pairs of iterations. 

%Starting all replicas from the volume-equilibrated con?guration described above, 100 iterations were conducted with 1 ps between exchange attempts to equilibrate the replicas to their respective temperatures. This equilibration run was followed by a production run of 500 iterations with 20 ps between exchange attempts, a total of 10 ns/replica. Solute configurations and potential energies from the production run were written to disk every 0.1 ps, while full-system restart files were recorded every 1 ps for the purpose of starting new simulations from these configurations, as described in section 3.4. 

%\subsection{State Decomposition (from MMS)}

%The slow degrees of freedom for terminally blocked alanine peptide (neglecting those involving solvent motion) can be captured by the two backbone torsion angles labeled $\phi$ and $\psi$ (see Figure 1) \cite{Bolhuis:2000p683,Ma:2005p819}. To this end, the potential of mean force at 302 K was computed from the parallel tempering data using the weighted histogram analysis method (WHAM) \cite{Kumar:1992p860,Chodera:2007p2346} and is shown in Figure 1. Six free energy basins are readily visible, and rectangular regions around these basins were chosen for the decomposition of all of configuration space into a set of six states. State definitions are listed in Table 1 and plotted as thick dividing lines in Figure 1. 
\end{widetext}
\end{document}