% PRL look and style (easy on the eyes)
\documentclass[aps,pre,twocolumn,showkeys,nofootinbib,superscriptaddress,linenumbers]{revtex4-1}
%\documentclass[pre]{revtex4-1}
% Double-spaced.
%\documentclass[aps,endfloats,prl,preprint,nofootinbib,superscriptaddress]{revtex4-1}

%\usepackage{palatino}
%\usepackage{newcent}
\usepackage{utopia} % font size about same
%\usepackage{bookman}
\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{{figures/}}

%\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}}}  % QED mark
\newcommand{\var}[1]{\mathrm{var} \, #1}
\newcommand{\cov}[2]{\mathrm{cov}\left( {#1}, {#2} \right)}

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

%% TITLE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{Supplementary Material for:\\
A robust approach to estimating rates from time-correlation functions}

\author{John D. Chodera}
 \email{jchodera@berkeley.edu}
 \affiliation{California Institute of Quantitative Biosciences (QB3), University of California, Berkeley, California 94720, USA}
 
\author{Phillip J. Elms}
 \email{elms@berkeley.edu}
 \affiliation{Biophysics Graduate Group, University of California, Berkeley, California 94720, USA}
 \affiliation{Jason L.~Choy Laboratory of Single Molecule Biophysics, University of California, Berkeley, CA 94720, USA}

\author{William C. Swope}
 \email{swope@us.ibm.com}
 \affiliation{IBM Almaden Research Center, San Jose, CA 95120}

\author{Jan-Hendrik Prinz}
 \email{jan.prinz@fu-berlin.de}
 \affiliation{DFG Research Center Matheon, Freie Universit\"{a}t Berlin, Arnimallee 6, 14195 Berlin, Germany}

\author{Susan Marqusee}
 \email{marqusee@berkeley.edu}
 \affiliation{Department of Molecular \& Cell Biology, University of California, Berkeley, CA 94720, USA}
 \affiliation{California Institute of Quantitative Biosciences (QB3), University of California, Berkeley, California 94720, USA}
 \affiliation{Jason L.~Choy Laboratory of Single Molecule Biophysics, University of California, Berkeley, CA 94720, USA}

\author{Carlos Bustamante}
 \email{carlos@alice.berkeley.edu}
 \affiliation{Department of Molecular \& Cell Biology, University of California, Berkeley, CA 94720, USA}
 \affiliation{California Institute of Quantitative Biosciences (QB3), University of California, Berkeley, California 94720, USA}
 \affiliation{Department of Physics, University of California, Berkeley, CA 94720, USA}
 \affiliation{Department of Chemistry, University of California, Berkeley, CA 94720, USA}
 \affiliation{Howard Hughes Medical Institute, University of California, Berkeley, CA 94720, USA}

\author{Frank No\'{e}}
 \email{frank.noe@fu-berlin.de}
 \affiliation{DFG Research Center Matheon, Freie Universit\"{a}t Berlin, Arnimallee 6, 14195 Berlin, Germany}

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

\date{\today}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MAKE FONT SMALLER TO MATCH JOURNAL STYLE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\small

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ABSTRACT/pacs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{abstract}
%\end{abstract}

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

\maketitle

\small % to match journal font size

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MODEL SYSTEM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Model system}

To illustrate how na\"{i}ve rate estimation by the crossing rate estimator $k_\mathrm{crossing}$ degrades as the reaction coordinate quality becomes poor, we consider a model two-state system, with states located at $x = 1$ and $x = 2$ and matrix of rate constants:
\begin{eqnarray}
\bfm{K} &\equiv& \left[ \begin{array}{rr}
- 0.052680 & 0.052680 \\
 0.052680 & - 0.052680
\end{array}\right]
\end{eqnarray}
such that the phenomenological rate constant $k = 0.105360$.
A discrete hopping trajectory of 50 000 observations was generated using an observation interval of $\Delta t = 1$, and Gaussian noise with zero mean and variance $\sigma^2$ was added to this discrete trajectory to mimic the effect of reaction coordinate degradation, with $\sigma$ chosen to be 0.1, 0.2, 0.3, or 0.4 (see Figure~\ref{figure:model-rate-estimation}). 
The larger $\sigma$, the poorer the quality of the reaction coordinate due to the difficulty in resolving which state the system is in near the transition region.

As can be seen in Figure~\ref{figure:model-rate-estimation}, the competing effects of a poor division between states by $x^\ddag$ and a finite observation interval $\Delta t$ result an \emph{underestimate} of the true rate by $k_\mathrm{crossing}$ for small $\sigma$ (dashed black line in middle column, top row), eventually becoming an \emph{overestimate} of the true rate for large $\sigma$ (dashed black line in middle column, bottom row).
The rate estimates from reactive flux $k_\mathrm{RF}(t)$ (dash-dotted black lines) are generally more close to the true rate and less sensitive to the location of the dividing surface, but the plateau region quickly drops off after $\tau = 2$ (right column).
The implied rate estimate $k_\mathrm{im}(t)$ (solid black line), however, produces a consistently good and relatively insensitive estimate of the true rate (which is shown in green).

\begin{figure*}[tbp]
\resizebox{\textwidth}{!}{\includegraphics{model-rate-estimation.pdf}} 
\caption{\label{figure:model-rate-estimation} {\bf Rate estimates for model two-state system.}
A two-state discrete state trajectory between two states centered on $x = 1$ and $x = 2$ was generated, on top of which was superimposed Gaussian noise with standard deviation $\sigma$ of 0.1 (top), 0.2 (upper middle), 0.3 (lower middle), or 0.4 (bottom).
\emph{Left:} The first 5 000 observations of the noise-contaminated observed trajectory $x(t)$ is shown on the left, with the optimal dividing surface $x^\ddag = 1.5$ drawn in red.  The shaded plot to the right of the trajectory depicts a histogram of the entire trajectory of 50 000 observations.
\emph{Middle:} Rate estimate from crossing rate $k_\mathrm{crossing}$ (red), reactive flux rate estimated $k_\mathrm{RF}(\tau)$ near plateau time of $\tau = 2$ (green), and implied rate $k_\mathrm{im}(\tau)$ evaluated at $\tau = 20$ (black), estimated as a function of dividing surface $x^\ddag$ choice the the vicinity of the transition region.
The true rate $k$ is shown as a dashed blue line.
\emph{Right:} Rate estimates from crossing rate $k_\mathrm{crossing}$ (red), reactive flux rate estimated $k_\mathrm{RF}(\tau)$ (green), and implied rate $k_\mathrm{im}(\tau)$ (black), evaluated at optimal dividing surface $x^\ddag = 1.5$.
The true rate $k$ is shown as a dashed blue line.
} 
\end{figure*}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ESTIMATION OF RATES AND STATISTICAL ERROR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Statistical error of implied rate estimate}
\label{appendix:statistical-error}

We first note that the implied rate estimate $k_{im}(t)$ given by Eq.~16 in the main text can be written as,
\begin{eqnarray}
k_{im}(t) &=& - t^{-1} \ln \frac{C_{AB}(t)}{\pi_A (1-\pi_A)} ,
\end{eqnarray}
where we have defined the correlation function $C_{AB}(t)$ as,
\begin{eqnarray}
C_{AB}(t) &\equiv& \expect{h_A(0) h_B(t)} .
\end{eqnarray}

To estimate the statistical error in the implied rate $k_\mathrm{im}(t)$, we first define two trajectory functionals:
\begin{eqnarray}
\mathcal{F}[X] &\equiv& h_A(x_0) \nonumber \\
\mathcal{G}[X] &\equiv& h_A(x_0) \, h_B(x_t) .
\end{eqnarray}
We define a timeseries $F_n$ and $G_n$ obtained from evaluating these functionals on sequential (potentially overlapping) segments $X_n$ of a much longer trajectory, or multiple independent short trajectories, depending on what kind of data has been collected,
\begin{eqnarray}
F_n &=& \mathcal{F}[X_n] \nonumber \\
G_n &=& \mathcal{G}[X_n] .
\end{eqnarray}
Estimates of their expectations are given by the sample means,
\begin{eqnarray}
\pi_A &\approx& \hat{F} = \frac{1}{N} \sum_{n=1}^N F_n\nonumber \\
C_{AB}(t) &\approx& \hat{G} = \frac{1}{N} \sum_{n=1}^N G_n \label{equation:sample-means}
\end{eqnarray}
We compute the rate constant for a given observation time $t$ (whose functional dependence we shall suppress) from an estimate of the second eigenvalue $\hat{\mu}_2$:
\begin{eqnarray}
\hat{k}_{im} &=& - t^{-1} \, \ln \hat{\mu}_2
\end{eqnarray}
The second eigenvalue is estimated from a modified form of Eq.~15 from the main text,
\begin{eqnarray}
\mu_2 &=& 1 - \frac{C_{AB}(t)}{\pi_A (1 - \pi_A)} \nonumber \\
\Leftrightarrow \:\: \hat{\mu}_2 &=& 1 - \frac{\hat{G}}{\hat{F} (1 - \hat{F})} .
\end{eqnarray}

The variance of the estimate $\hat{k}_{im}(t)$ can be estimated by simple first-order Taylor series expansion propagation of error,
\begin{eqnarray}
\delta^2 \hat{k}_{im} &=& \left[ \frac{\partial \hat{k}_{im}}{\partial \hat{\mu}_2} \right]^2 \, \delta^2 \hat{\mu}_2 = t^{-2} \frac{\delta^2 \hat{\mu}_2}{\hat{\mu}_2^2}
\end{eqnarray}
We apply the first-order Taylor series expansion propagation of error to compute the uncertainty in $\delta^2 \hat{\mu}_2$:
\begin{eqnarray}
\delta^2 \hat{\mu}_2 &=& \left[ \frac{\partial \hat{\mu}_2}{\partial \hat{F}} \right]^2 \delta^2 \hat{F} + \left[ \frac{\partial \hat{\mu}_2}{\partial \hat{G}} \right]^2 \delta^2 \hat{G} + 2 \left[ \frac{\partial \hat{\mu}_2}{\partial \hat{F}} \right] \left[ \frac{\partial \hat{\mu}_2}{\partial \hat{G}} \right] \delta \hat{F} \delta \hat{G} \nonumber \\
\end{eqnarray}
where the required derivatives are given by
\begin{eqnarray}
\frac{\partial \hat{\mu}_2}{\partial \hat{F}} = \frac{\hat{G} (1 - 2 \hat{F})}{\hat{F}^2 (1 - \hat{F})^2} \:;\: \frac{\partial \hat{\mu}_2}{\partial \hat{G}} = \frac{1}{\hat{F} (1 - \hat{F})} 
\end{eqnarray}
To estimate $\delta^2 \hat{\mu}_2$, we must estimate the variance and covariance of the estimators $\hat{F}$ and $\hat{G}$:
\begin{equation}
\delta^2 \hat{F} = \frac{\var{F_n}}{N/g} \:;\: \delta^2 \hat{G} = \frac{\var{G_n}}{N/g} \:;\: \delta \hat{F} \delta \hat{G} = \frac{\cov{F_n}{G_n}}{N/g}
\end{equation}
where the sample covariances are used to estimate $\var{A_n}$, $\var{B_n}$, and $\cov{A_n}{B_n}$, and $g$ is the maximum statistical inefficiency for the timeseries $F_n$ and $G_n$.
If the $X_n$ denote independent short trajectories (such as obtained by a computer simulation from uncorrelated initial starting configurations $x_0$), then $g = 1$; otherwise, the statistical inefficiency $g$ can be estimated via standard means (see Section~5.2 of \textcite{chodera:jctc:2007}).

Practically, we take advantage of time-reversibility of dynamics, and use a slightly modified set of trajectory functionals that yield the same expectation but average over more snapshots from the trajectory in the case that the trajectory segments are of length $T > t$:
\begin{eqnarray}
\mathcal{G}[X] &=& \frac{1}{T-t} \sum_{t_0 = 0}^{T-\tau} \frac{1}{2} \left[ h_A(x_0) \, h_B(x_{t_0+t})  + h_B(x_{t_0}) \, h_A(x_{t_0+t}) \right] \label{equation:trajectory-functionals-for-rate} \nonumber \\
\mathcal{F}[X] &=& \frac{1}{T-t} \sum_{t_0 = 0}^{T-t} \frac{1}{2} \left[ h_A(x_0)  + h_A(x_{t_0+t}) \right] 
\end{eqnarray}
This appropriately accounts for the fact that all time origins produce equally valid estimates and, for systems with multiple conformational states, ensures satisfaction of detailed balance.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SCRIPTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Matlab code}

All Matlab scripts used to compute rate estimates in this paper, as well as all datasets and plotting scripts, are available online at \url{http://simtk.org/home/rate-theory/}.

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

\bibliographystyle{apsrev4-1} 
\bibliography{rate-theory}

\end{document}