% PRL look and style (easy on the eyes)
\documentclass[twocolumn,english,aps,pre,nofootinbib,superscriptaddress]{revtex4}
% Two-column style (for submission/review/editing)
%\documentclass[aps,prl,preprint,nofootinbib,superscriptaddress]{revtex4}
\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
\usepackage{amsmath}
\usepackage{amssymb}

\makeatletter

\usepackage{palatino}
%\usepackage[mathbf,mathcal]{euler}
%\usepackage{citesort}
\usepackage{graphicx}
\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


\makeatother

\usepackage{babel}


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



%% TITLE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{Optimal use of data in parallel tempering simulations for the construction of kinetic models of peptide dynamics}


\author{Jan-Hendrik Prinz$^{1}$}
 \email{jan-hendrik.prinz@iwr.uni-heidelberg.de}
 \affiliation{IWR, University of Heidelberg, INF 368, 69120 Heidelberg, Germany}
 \thanks{$^{1}$\textbf{ }Both authors have equally contributed to this work}

\author{John D. Chodera$^{1}$}
 \thanks{Corresponding author}
 \email{jchodera@stanford.edu}
 \affiliation{Department of Chemistry, Stanford University, Stanford, CA 94305}


\author{Vijay S. Pande}
 \email{pande@stanford.edu}
 \affiliation{Department of Chemistry, Stanford University, Stanford, CA 94305}

\author{Jeremy C. Smith}
 \email{smithjc@ornl.gov}
 \affiliation{ORNL, Oak Ridge, TN 12345}

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

\date{\today}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ABSTRACT/pacs 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}
	- Recently, others have demonstrated how the short physical trajectories generated in parallel tempering simulations between exchanges can be used to construct Markov models to describe peptide dynamics at each simulated temperature. 
	
	- While able to describe the temperature-dependent kinetics, this approach does not make optimal use of all available data, instead parameterizing the model at a given temperature with only data from that temperature. 
	
	- We propose a strategy where, with a simple modification of the parallel tempering protocol, the harvested trajectories can be reweighted, allowing data from all temperatures to contribute to the estimated kinetic model without requiring the assumption of a particular kinetic model (like Arrhenius) beyond that of classical statistical mechanics. 
	
	- This method not only reduces the statistical uncertainty in the kinetic model, providing estimates of transition rates even for transitions not observed at the temperature of interest, it allows the kinetic model to be estimated at temperatures other than those simulated from. 
	
	- To demonstrate the method in a system, we apply the method to the solvated terminally-blocked alanine peptide.
\end{abstract}

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

\maketitle

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

Conformational transitions are critical to the function of proteins and nucleic acids. These transitions span large ranges of lengthscales, timescales and complexity, and include ligand binding \citet{Nienhaus_Nature00}, complex conformational rearrangements between native protein substates \citet{FischerSmith_PNAS102_6873,NoeKrachtusSmithFischer_JCTC}, and folding \citet{JaegerEtAl_PNAS06,KobitskiNienhaus_NucleicAcidsRes07}. Understanding these processes is challenging, as they often involve various pathways \emph{via} many intermediate conformations. 

A natural approach towards modeling the kinetics of molecules is by partitioning the state space into discrete states \citet{Wales,NoeHorenkeSchutteSmith_JCP07,ChoderaEtAl_JCP07}. In particular, biomolecular function often depends on the ability to undergo transitions between long-lived, or {}``metastable'' states \citet{Frauenfelder_Science254_1598}, which serve as a useful definition of states in a kinetic model \citet{NoeHorenkeSchutteSmith_JCP07,ChoderaEtAl_JCP07}.

The switching process between conformational substates is often described with the memoryless Master equation:

\begin{equation}
	\frac{dp(t)}{dt}=p(t)K.
	\label{eq_Master-equation}
\end{equation}

with $p(t)\in\mathbb{R}^{1\times m}$ containing the probability to find the system in each of its $m$ states at time $t$. $K\in\mathbb{R}^{m\times m}$ is a rate matrix with $K_{ij}$ being the transition rate from state $i$ to state $j$ and the diagonal elements are given by $K_{ii}=-\sum_{j\neq i}K_{ij}$ in order to ensure mass conservation. Alternatively, the system dynamics can be described by a discrete-time Markov process using the transition matrix, $T(\tau)\in[0,1]^{m\times m}$, whose entries, $T_{ij}$, provide the probability of the system to be found in state $j$ at time $t+\tau$ given that it was in state $i$ at time $t$. $T(\tau)$ is row-stochastic, \emph{i.e.}: $\sum_{j}T_{ij}=1$ for all $i$.

The time-discrete analog to Eq. (\ref{eq_Master-equation}) is:

\begin{equation}
	p(k\tau)=p(0)T^{k}(\tau).
	\label{eq_Chapman-Kolmogorov}
\end{equation}

Eq. (\ref{eq_Master-equation}) and (\ref{eq_Chapman-Kolmogorov}) provide equivalent results at discrete times $t=k\tau$, $k\in\mathbb{N}_{0}$ and are related by $T(\tau)=\exp(\tau K)$ \citet{vanKampen}. Here, the transition matrix formulation will be used. 

\begin{itemize}
	\item In order for the dynamics to be markovian at lagtime $\tau$ states need to be defined such, that their internal equilibraiation time is significantly shorter than the lagtime $\tau$. This is not trivial to ensure for complex molecules, however this is not the issue of the current article and is discussed elsewhere. Cite: (Hummer), JCP(John), JCP(Frank, Illia, s.o)
	\item It has been shown, that using appropiate definitions of the discrete states the essential dynamics of biomolecules is well-described by such a markov model. Cite: MMS(John), JCP07(Frank, Illia, ...), JPCB?08 (Hummer).	
	\item The estimation of the transition probabilities is hampered by the sampling problem, due to which transition between meta-stable states are rare events. It would be desirable to use dynamics of this same system at some higher temperature where these transitions occur more easily / frequently
	\item Parallel tempering simulations and Replica Exchange Molecular Dynamics(Cite: PT...), in short PT simulations, have become popular means to explore the state space at low temperatures using a set of coupled simulations at a range of temperatures. The PT simulations at sifferent temperature are coupled by exchanges constructed in a way that they converge toward a canonical distribution at each temperature. PT simulations have been widely used to explore conformations in biomolecules, that are difficult to access via direct MD.
	\item Recently, an approach to estimate interstate transition probabilities from the short physical trajectories generated from PT trajectories between exchanges has been proposed. (Cite: PRE (Hummer)). However, this approach produces a separate estimate of the transition matrix at each temperature without attempting to combine the data at all available temperatures.
	\item We propose an approach to obtain an optimal estimation of the temperature dependend transition matrix combing data from all PT temperatures. The estimated tranition matrix is statistically optimal in a maxmimum likelihood sense and does not require any model assumptions such as the Arrhenius Law. the estimation procedure is based on \citet{shirts-chodera:jcp:2008:mbar}. A slight modification of the PT protocol is necessary.
	\item As a result of the optimization procedure one obtains the estimate of the transition matrix at any temperature desired. Even at temperatures, that have previously not been sampled from. Additionally an estimate for the error of the tranision matrix elements is obtained.
	\item The procedure is illustrated on a solvated terminally-blocked alanine dipeptide simulated with PT at temperature in the range 273K upto 600K.
	\item Outline of this manuscript
\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THEORY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theory}
\label{section:theory}

\subsection{Discrete-state Markovian kinetic models}

Let $s_{i},i\in\left\{ 1,\ldots,n\right\} $ be a complete ($\bigcup_{i=1}^{n}s_{i}=\Omega$) and non-overlapping ($s_{i}\cap s_{j}=\emptyset,\forall i\ne j\;\big|\; i,j\in\{1,\ldots,n\}$ ) partition of state space $\Omega$. These $s_{i}$correspond to the state between which interstaet transition probabilites will be computed.

Define the characteristic function $\chi$ as

$\chi_{i}(x)=\begin{cases}
		1 & \text{if }x\in s_{i}\\
		0 & \text{else}
	\end{cases}.$

{\color{red} 

The transition probability from state $i$ to state $j$ at lagtime $\tau$ and inverse temperature $\beta$ is then defined as 
\begin{eqnarray}
	T_{ij}(\tau;\beta) & = & \frac{C_{ij}(\tau;\beta)}{\sum\limits _{k}C_{ik}(\tau;\beta)}
	\label{equation:transition-matrix-definition}
\end{eqnarray}

defining the row-stochastic transition matrix $T(\tau,\beta)\in\mathbb{R}^{n\times n}$ Cite: (Bill). The state-to-state time-correlation function $C_{ij}(\tau;\beta)$ is given by 
\begin{eqnarray}
	C_{ij}(\tau;\beta) & \equiv & \expect{\chi_{i}(0)\,\chi_{j}(\tau)}_{\beta}
\end{eqnarray}  
and is symmetric --- $C_{ij}(\tau;\beta)=C_{ji}(\tau;\beta)$ --- for systems obeying detailed balance. }

\subsection{Transition probabilities by reweighting}

\subsubsection{Dynamical reweighting}

Suppose we have a set of $N_{k}$ Hamiltonian trajectories $\z_{kn}(t)$, $n=1,\ldots,N_{k}$, $t\in[0,T]$, where 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}~\citet{shirts-chodera:jcp:2008:mbar}, we can estimate a correlation function $C_{AB}(\tau;\beta)$ using this set of trajectories as 
\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))$. The normalized trajectory weights $w_{kn}(\beta)$ are specified in terms of unnormalized weights $\tilde{w}_{kn}(\beta)$ 
\begin{eqnarray}
	\tilde{w}_{kn}(\beta) & = 		& \left[\sum_{k'=1}^{K}N_{k'}\,\exp[f_{k'}-(\beta-\beta_{k'})\, E_{kn}]\right]^{-1}
\end{eqnarray}

such that $w_{kn}(\beta)=\tilde{w}_{kn}(\beta)/\sum\limits _{k=1}^{K}\sum\limits _{n=1}^{N_{k}}\tilde{w}_{kn}(\beta)$, with $E_{kn}\equiv H(\z_{kn}(0))$ denoting the total energy of the trajectory.

The dimensionless free energies $f_{i}=-\ln Z_{i}+c$ are determined by solution of a set of self-consistent equations 
\begin{eqnarray}
	e^{-f_{i}} & = & \sum_{k=1}^{K}\sum_{n=1}^{N_{k}}\tilde{w}_{kn}(\beta)
\end{eqnarray}
which can be done efficiently by a number of means~\citet{shirts-chodera:jcp:2008:mbar}.

%The statistical uncertainty in the estimate of $C_{AB}(\tau;\beta)$ can also be estimated by forming the matrix $\bfm{\Theta}$:
%\begin{eqnarray}
%\Theta_{ij} &=&
%\end{eqnarray}

\subsubsection{Modified parallel tempering protocol}

By using a modified parallel tempering protocol, we can generate a set of Hamiltonian trajectory segments $\z_{kn}(t)$ of uniform length $T\ge\tau$ whose initial phase space points $\z_{kn}(0)$ are distributed from the canonical (NVT) ensemble at corresponding temperatures $\beta_{1},\ldots,\beta_{K}$.
We start assuming that some process was used to generate 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 H(\z_{k0}(0))}
\end{eqnarray}
This could be done by means of a standard PT 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$, we propagate Hamilton's equations of motion using a symplectic integrator to generate trajectories of $\z_{kn}(t)$ of length $T$. Finally, we propose exchanges between the final configurations $\z_{in}(\tau)$ and $\z_{jn}(\tau)$ of temperatures $\beta_{i}$ and $\beta_{j}$, accepting or rejecting the exchange with the Metropolis-like probability depending 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\} 
\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 as $\z_{k(n+1)}(0)$, from which the next iteration can be carried out.

\subsubsection{Transition probabilities}

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_{knij}(\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_{knij}(\tau) 	& \equiv 	& \frac{1}{2(T-\tau)}\int_{0}^{T-\tau}dt_{0}\,[\chi_{i}(\z_{kn}(t_{0}))\,\chi_{j}(\z_{kn}(t_{0}+\tau))\nonumber \\
					&  			& \mbox{}+\chi_{j}(\z_{kn}(t_{0}))\chi_{i}(\z_{kn}(t_{0}+\tau))]
	\label{equation:definition-A-matrix}
\end{eqnarray}
The row-stochastic transition matrix estimate $\hat{\bfm{T}}(\tau;\beta)$ is then estimated from Eq.~\ref{equation:transition-matrix-definition}
\begin{eqnarray}
	\hat{T}_{ji}(\tau;\beta) & = & \frac{\hat{C}_{ij}(\tau;\beta)}{\sum\limits _{k}\hat{C}_{ik}(\tau;\beta)}
\end{eqnarray}
Because $\hat{\bfm{C}}(\tau;\beta)$ is symmetric, $\hat{\bfm{T}}(\tau;\beta)$ will be reversible.

\subsection{Bayesian estimate of transition probabilities from a single temperature}
\label{subsection:transitionmatrix-sampling}

\begin{itemize}
	\item To compare the estimates of the present method 
	\item To analyze, how the combiantion of data from all temperatures reduces the uncertainy in the estimation of the transition matrix we compare with a method to compute errors at a single temperature. The method is described in detail in Cite: JCP08(Frank) and is briefly outlined here.
	\item Let $C_{ij}(\tau,\beta)$ be the transition counts from all trajectories at inverse temperature $\beta$.Then the posterior probability of transition matrices given this observation is
\end{itemize}

\begin{equation}
	p(T|C)\propto p(T)p(C|T)=p(T)\prod_{i,j\in S}T_{ij}^{c_{ij}},
	\label{equation:Posterior-relation}
\end{equation}

We use a uniform prior distribution $p(T)\equiv const$ and thus $p(T|C)\propto p(C|T)$ Furthermore we restrict ourselves to transition matrices, which fulfill detailed balance, i.e. are reversible with respect to the stationary distribution $\bfv{\pi}$:

\begin{equation}
	\pi_{i}T_{ij}=\pi_{j}T_{ji}.
\end{equation}

Distribution (\ref{equation:Posterior-relation}) was sampled using a Markov chain Monte-Carlo procedure described in Cite: JCP08(Frank).


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

\begin{itemize}
	\item Description of system used in the test: solvated terminally-blocked alanine peptide, relevant details from the protocol (Copy from MMS (John))
	\item PT protocol that was used. Integrator, energy conservation - validated, fluctuations, drift below ?, temperatures, burn-in phase / equilibration 1ns. to ensure, that all initial configurations are drawn from equilibrium at their respective temperatures. 40 temperatures exponentially spaced between 273K and 600K, 500 trajectories each with a length of 20ps saved every 100fs
	\item State definition $\chi_{i}$: 6 states have been manually defined in the $\phi,\psi$ space of dihedrals. Using this state-definition the interstate dynamics have been shown to be approximately Markovian at a lagtim of $\tau=6ps$.\citet{chodera:mms:2006}, 2004a or b(Swope) Picture -> Energy Landscape for Peptide in $\phi,\psi$ space. The same state definition was used across all temperatures.
	\item Each contiguous trajectory snippet $n$ at temperature $k$ the corresponding symmetrized(reversible) state-to-state time-correlation function $A_{knij}$ was computed using Eq. (\ref{equation:definition-A-matrix})
	\item Determination of the weights $\tilde{w}_{kn}(\beta)$, involving estimation of {}``freeenergies'' $f_{k}$ . From this computing transition probabilites $\hat{T}_{ij}$and related uncertainties(should be contained in the theory/methods section, if too long move parts to appendix)
	\item Comparison of transition probabilities as a function of temperature. (Picture) Description of Plots {}``Blue lines show the estimates Transition Probabilities using the present method, the result. Error bars indicate 2$\sigma$ confidence intervals. For comparison the black lines show the maximum-likelihood transition probabilities computed from single temperature data. It is apparent, that the reweighting methods provides useful and bounded estimated across the whole temperature range even when the corresponding transition was not observed at that temperature. The counting method without reweighting is very noisy and not useful in situations at temperature with no transitions counted.
	\item Red lines show Transition Probabilities using reversible TM Sampling conducted seperately at each individual temperature using the method described in section \ref{subsection:transitionmatrix-sampling}. 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 contraint. Apparently the reweighting estimate that uses the combined data from all temperatures has smaller uncertainties than the estimators that work on individual temperatures.
	\item For transitions, that are not sampled are certain temperature ranges, the maximum-likelihood estimate obtained with our reweigthing method is approximately zero(s. fig ... panel ...) 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 ...)
	\item To validate our result. Further simulations were conducted at $T=302K$. 10000 trajectories shoot from the equilibrium of each state with length 10ps. Resulting in a single transition at a lagtime $\tau=6ps$. The results of these simulations are plotted in figure ... as the black line. The blue curve is the reweighting estimate using combined date from all temperatures of the short PT simulation and the gray histogram the sampled distribution using only single temperatures from the PT simulation. In many cases the combined temperature estimate is closer to the reference distribution than the sampled distribution.	
	\item Comparison of temperature dependence of implied timescales -Description. s. fig. ...
	\item Some figure showing how much different temperatures contribute to the estimates of particular transitions - Description. s fig. 
	\item In Extrapolation figure? Is this reasonable at all ?
\end{itemize}



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

\begin{itemize}
	\item Even at temperatures at which no transitions are observed, can still estimate transition probabilities. 
	\item Allows us to interpolate at temperatures that haven't been sampled; smoothly describe temperature-dependent behavior (fluctuations much smaller than Bayesian method) 
	\item We could even differentiate the estimates of transition probabilities if we wanted because trajectory weights are differentiable 
	\item Same is true for non-dynamic properties (expectations, stationary distributions) 
	\item Some transitions appear to be thermally activated processes (increasing probability with energy or temperature), others are rather flat 
	\item What about entropically-dominated barriers? In this case, heating would reduce rate, but because we can use lower temperatures, this can help too. 
	\item While temperature range at which other temperatures contribute significantly to weighting diminishes as system size increases, the requirement that temperatures in parallel tempering be spaced to get good exchange means that several replicas are always guaranteed to contribute 
	\item Any observations about metastable state structure or grouping as temperature changes? (PCCA+ as a function of temperature?) 
	\item What about the future? Some combination of Bayesian and reweighting (similar to Ron Levy's T-WHAM) may provide the best of both types of estimators in giving good uncertainties at the expense of introducing some bias from introduction of histograms in energy 
	\item Transition path sampling could help collect more data for poorly sampled transitions -- comment on \char`\"{}sweet spot\char`\"{} for energies of these transitions? 
	\item Limitations of requiring we use Hamiltonian trajectories -- could it work for longer trajectories or non-Hamiltonian dynamics?
\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ACKNOWLEDGMENTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Acknowledgments}
\label{section:acknowledgments}

JDC and VSP gratefully acknowledge support from an NSF grant for Cyberinfrastructure (NSF CHE-0535616).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BIBLIOGRAPHY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{1_Users_jan-hendrikprinz_Studium_Projekte_dynam___kov-model-from-parallel-tempering_lyx_prsty}
\bibliography{0_Users_jan-hendrikprinz_Studium_Projekte_dynam___ng_lyx_markov-model-from-parallel-tempering}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% APPENDICES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\appendix


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% APPENDIX: PROOF THAT MODIFIED PARALLEL TEMPERING PROTOCOL GENERATES CORRECT DISTRIBUTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Proof that modified parallel tempering protocol generates canonical distribution}
\label{appendix:parallel-tempering-proof}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% APPENDIX: PROPAGATION OF ERROR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Estimation of uncertainties in transition probabilities}
\label{appendix:propagation-of-error}



\end{document}
