\documentclass{revtex4-1}

\usepackage{palatino}
\usepackage{amsmath}
\usepackage{amssymb}

\usepackage{graphicx}
\usepackage{dcolumn}
\usepackage{boxedminipage}
\usepackage{verbatim}
\usepackage[colorlinks=true,citecolor=blue,linkcolor=blue]{hyperref}

% The figures are in a figures/ subdirectory.
\graphicspath{{figures/}}

% 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{\bfc}{\bfm{c}}
\newcommand{\hatf}{\hat{f}}

\newcommand{\bTheta}{\bfm{\Lambda}}
\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}_{#1}\left[ #2\right]}
\newcommand{\T}{\mathrm{T}}                                % T used in matrix transpose

% Commands to force LaTeX to put two figures on a page:
\renewcommand{\textfraction}{0.05}
\renewcommand{\topfraction}{0.95}
\renewcommand{\bottomfraction}{0.95}
\renewcommand{\floatpagefraction}{0.35}
\setcounter{totalnumber}{5}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DOCUMENT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\begin{center}
{\bf  Responses to reviewers:}
\end{center}

We thank the reviewers for their kind words about this paper and the
suggestions they had for clarifications.  We have revised our
manuscript accordingly, setting the revisions in red text.  Below, we
address the referee comments individually; original referee comments
are set in blue text.  We have also fixed a small error in the
``restricted range'' variant of the algorithm that we caught while the
manuscript was in review and made a number of small stylistic changes
for improved clarity.

{\color{blue} \bf Reviewer \#1:}
 
{\color{blue} (1) The authors have sidestepped the issue that limits the utility of
all expanded ensemble methods, namely that in order to be effective, a
sampleable ensemble needs to be identified. This in no way detracts
from the contribution of the paper, which concerns improving the
mixing of coupled Markov chains, but perhaps deserves at least a
mention in the introduction.}

We are not entirely clear what the referee means by ``a sampleable
ensemble needs to be identified,'' but we assume this has to do with
the fact that it may not always be possible or easy to identify a set
of states with disjoint state variables that can be treated by the
expanded ensemble method.  Certainly, whenever there is a group of
distributions differing by some parameter affecting the Hamiltonian
(such as temperature, or the state of the alchemical Hamiltonian),
then we will be able to construct an expanded ensemble that we can set
and find weights for (another difficult challenge that we do not
directly address in this paper).  We are assuming that the readers
will generally be able to define such an expanded ensemble for their
system (and we also cite literature on various schemes for doing this
optimally), and focus on how they can accelerate the sampling given
such an expanded ensemble.

%For any specified set of log weights $\{g_k\}$, it is in principle possible to sample the resulting ensemble described by the corresponding probability density $\pi(x,k) \propto \exp[-u_k(x) + g_k]$, though it may sometimes be challenging to do so.
%Our manuscript presents a simple way to reduce the effort required to generate uncorrelated samples from this distribution by making simple modifications to one step of the sampling procedure.

%The referee \emph{may} be referring to the process by which a set of log weights $g_k$ are determined.
%We note that, though the most common approach is to attempt to select the weights to ensure approximately equal sampling of each state $k$, such that $g_k = f_k \equiv - \ln \int e^{-u_k(x)}$, this is not always the case.
%Other choices of weights may be desirable, depending on the ensemble one wishes to sample.
%For example, if the index $k$ denotes a ligand species for a protein receptor, selecting $g_k$ to be an approximation to the solvation free energy of the ligand in the absence of the receptor will bias the ensemble toward ligand species that bind more tightly to the receptor, which may be desirable.
%Even if one wishes to use the so-called ``perfect weights'' $g_k = f_k$, there are so many different strategies for choosing these weights adaptively or through pilot simulations that it would be far beyond the scope of this work to review them in our manuscript, let alone evaluate them.
%We expect our improvements to carry over to many of these adaptive schemes quite generally.

%A discussion of this issue, as well as numerous references, is already given in Section II.C of the manuscript, and our justification for fixing the weights to ``near perfect'' values for purposes of simplifying our analysis was set out in Section V.

{\color{blue} (2) Is there any way to get a rough idea of how the improvements in
statistical efficiency will depend on system size and complexity?}

Sampling efficiency will always depend on the system and specifics of the simulation (such as the number of thermodynamic states $K$ and their definition), so it can be difficult to answer this question with generality.
In the case that configurational sampling is rapid compared to the time between exchange attempts, however, the efficiency gain will not depend on the system size or complexity, but rather the probability that exchange attempts between non-neighboring states will be accepted.  
%The key to the approaches discussed in this study is that it increases the movement through the state dimensions relative to nearest neighbor exchange.  
%The size of the system is not the determining factor.   
%If the movement in configurations space is significantly slower than that in state space, then the advantage of independence sampling will be relatively low. 
%In the end, this is not a cure all---it merely allows for faster sampling in coordinates that only used nearest neighbor sampling.  

{\color{blue} After all, ultimately one would like to use such methods for
calculations involving tens of thousands of atoms, to study protein
conformational change or ligand binding.}

There are a number of systems of interest or varying size and complexity, and given red the current length and general nature of our conclusions, it probably beyond the scope of the current manuscript to explore too many specific cases.  Protein conformation changes are certainly of particular interest to us.  We are currently investigating how much independence sampling techniques can improve the structural equilibration of small peptides using parallel tempering beyond what is obtained with standard nearest neighbor methods, but believe that the systems we examine in the current manuscript are better for introducing the ideas to the field.

{\color{blue} \bf Reviewer \#2:}

{\color{blue} The authors describe some simple modifications for accelerating the mixing of
thermodynamic states during replica exchange and expanded ensemble
simulations. They show that swapping thermodynamic states at random leads to a
better sampling compared to standard swapping approaches such as nearest
neighbor swaps in replica exchange simulations. Model-systems on which they
demonstrate this include a large Lennard-Jones sphere and an alanine
dipeptide. The manuscript is well written and both sections reflect a deep
understanding of the problem.}
 
{\color{blue} (1) Several alternative modifications to improve mixing are
  suggested. Even if the authors describe in details two methods easy
  to implement (Independence Sampling and Metropolized Indep.), they
  also show the equivalence with a third method that does not require
  any modification at all to the common algorithm ('nearest neighbor'
  repeated many times in a row). The simplest (third) proposed
  modification will probably be more appealing to users of replica
  exchange for example. In addition the reader may not understand, still
  after the method section, why/how Independence Sampling could ever
  increase the acceptance rate for exchanges of thermodynamic states
  compared to consecutive passes of the neighbor exchange algorithm
  since the latter was set up (4-10, 74) for a maximum acceptance
  rate. On larger systems, users of replica exchange will know that
  random swaps (exchange attempts between non-neighbors) lead to
  unacceptably small acceptance rate. The authors should investigate or
  at least discuss this issue.}

This is an important point to make clearly, and we thank the reviewer for helping us recognize that the original manuscript was not entirely transparent on this topic.

We are a bit confused by the specific way the reviewer has stated the
question; specifically, the statement that ``Independence Sampling
could ever increase the acceptance rate for exchanges of thermodynamic
states compared to consecutive passes of the neighbor exchange
algorithm''.  It's not entirely clear what ``consecutive passes''
means in this context.  We do claim that one can reach independence
sampling in the limit of a large number of sweeps of neighbor exchange
without any configuration sampling between them.  so if ``consecutive
passes'' means ``many sweeps of neighbor exchange'', then the methods
are indeed essentially identical.

The choice of spacing that maximizes neighbor exchange is still ideal
for independence sampling.  The key to independence sampling is that a
configuration that has nonnegligible probability in two to three
states more distant has the chance to move two or three states before
configurational updating occurs, as both types of sampling occur;
nearest neighbor and longer exchanges.

In order to explain this better, we have added to the paper the
sentence ``We note that the improvement in efficiency over neighbor
exchange algorithms is because at each state, there is a nonzero
probability of moving to all other states in proportion to the
Boltzmann weight difference, thereby including exchange with more
remote states without significantly downweighting exchange with
neighboring states.''


{\color{blue} (2) The fundamental problem with the approach proposed in this
  manuscript is that it has not been tested on complex systems with
  real application in biology. }

While the performance of these algorithmic modifications on large,
complex biomolecular systems is certainly of great interest and
importance, we feel that evaluation on these systems would be beyond
the scope of this first paper introducing the theory.

In our manuscript, we have chosen relatively simple systems to
characterize the behavior of our proposed algorithmic modifications:
specifically, a one-dimensional illustrative example and three
molecular systems (insertion of small and large Lennard-Jones spheres
in explicit solvent, parallel tempering of the alanine dipeptide in
implicit solvent, and 2D replica-exchange umbrella sampling of the
alanine dipeptide in implicit solvent).  We have chosen systems of limited
complexity in order to precisely characterize the effect of our
algorithmic modifications in reducing the overall correlation times
for generating fully independent samples, which required we run the
simulations sufficiently long to be able to show the difference in
correlation times with a high statistical significance.

We are in the process of carrying out further simulations to examine
whether and how methods such as independence sampling can accelerate
sampling for peptide folding.

{\color{blue} The authors mention that typical correlation time for water
  orientation is a few picoseconds, but this becomes much longer for
  non-trivial coordinates that result from collective moves during
  molecular dynamics. The authors mention that these collective moves
  might be disrupted by too frequent changes in state space. It would
  be great to know if/when this is the case.  Systematic jumps from
  one thermodynamic state to the other should prevent some key
  biomolecular processes to ever happen if (i) they have a slow
  relaxation time and (ii) they occur in a narrow range of values of
  the thermodynamic variables in question. The authors should
  investigate or carefully discuss this issue.}

The issues of systematic jumps disrupting collective motions is
essentially the same for any expanded ensemble or replica exchange
method that allows discrete transitions between potential energy
surfaces.  In this sentence, we were merely pointing out that this
could potentially be an issue with any expanded ensemble or replica
exchange method, and our current method should neither solve or
exacerbate this problem. In thinking through this problem more
completely, it is clear that things are a bit more
complicated---changes in the underlying potential energy surface can
significantly change the collective motions of the system, but there
is no simple proof that says that the new collective motions will be
either slower or faster.  We have therefore removed this section of the text.
 
{\color{blue} (3) "In the cases that proposal probabilities are based on past history
however, the algorithm will not preserve the equilibrium distribution" (p.4,
column 2):  Later (p.5) the Metropolized Indep. scheme is described as a particular case
of Independence Sampling where the current state is removed from the subset of
allowed states for the next move. What then prevents from removing some
previously visited states in addition to the current one? }

A transition probability that depends only on the current state (as
the Metropolized independence sampler does) is history independent and
Markovian---it depends only on the current state, not prior history.
Therefore, excluding the current state is still only dependent on the
current state.

{\color{blue} Some recent work (Bussi 2006, Roitberg 2007) indicates that
history-based generalized ensemble algorithms can be in equilibrium
under particular setting (e.g.~infinitesimally small bias combined
with stochastic replica exchanges, Bussi 2006). It might be better to
remove this statement or else develop better the conditions of
ergodicity under which Metropolized Indep. works.}
 
{\color{blue} \it Bussi G, Gervasio FL, Laio A, Parrinello M (2006) Free-energy
landscape for hairpin folding from combined parallel tempering and
metadynamics. J.  Am. Chem. Soc. 128: 13435-41.}

The history-dependent bias we discuss is exactly the "error in the
free-energy estimation is related to the diffusion coefficient in the
CV's dynamics" discussed by Bussi et al.  This error can be minimized
by choosing Gaussians such that the sampling in the CV is slower than
all other collective variables, but it is still there.  The Bussi 2006
paper from PRL discusses how to minimize this error, but it still
exists and has to be managed. We have now included this PRL citation.
 
{\color{blue} \it Roitberg AE, Okur A, Simmerling C (2007) Coupling of replica
  exchange simulations to a non-Boltzmann structure
  reservoir. J. Phys. Chem. 111: 2415-18. }

The history dependence of the reservoir is a bit
complicated. Reservoir methods are a very interesting new development
and the history dependence of the reservoir can be overcome, but as
with metadynamics, overcoming the errors due to history dependence is
also challenging.  These challenges are not really the sorts of things
we are addressing in the paper, and it would probably be best not to
confuse the reader by spending too much time discussing them.  In the
limit the reservoir is complete, sampling from the reservoir is a form
of indepenence sampling in configuration space (rather than state
space) with the potential to reduce correlation times, but it's not
clear that is what the reviewer is referring to.

%{\color{red}[JDC: I don't understand what you are trying to say here, or what the reviewer is trying to say.  Reservoir simulations do generally have a history-dependence, but they're supposed to be approximations to a different form of independence sampling in which uncorrelated samples from the reservoir seed the simulation and try to reduce correlation times.  I don't believe there has been any sort of detailed mathematical analysis of what kinds of problems the approximations most people actually use by using a finite or static reservoir introduce.]}

To attempt to address the reviewers comments, we have added a
reference to reservoir methods when discussing independence sampling
in configurations, and we have changed the sentence to read "In the
cases that proposal probabilities are based on past history however,
the algorithm will not preserve the equilibrium distribution, though
in some cases the deviations caused by the history dependence can be
mitigated with proper choices of parameters (Cite Bussi, G.; Laio, A.;
Parrinello, M. Phys. Rev. Lett. 2006, 96, 090601)

{\color{blue} (4) Typos:}

We thank the reviewer for taking the time to find and point out these typos.
They have been corrected in our revised manuscript.

%X - p.3, column 2&#x00A0;: closing bracket
% after&#x00A0;&#x00AB;&#x00A0;walker&#x00A0;&#x00BB; not opened.
%X - p.4, column 2&#x00A0;: (...) to sample in k or x depending according to
%(...). 
% - p.4, column 2&#x00A0;: closing bracket after ref. 71 not opened.
%X - p.5, column 2&#x00A0;: This scheme has the surprisingly property that,
%(...). 
%X - p.6, column 1&#x00A0;: (...) even something simple as application of (...).
%X - p.7, column&#x00A0;2 : (...) an expanded ensembles (...).
%X - p.7, column 2&#x00A0;: &#x00AB;&#x00A0;Nij is the number of times
%(...)&#x00A0;&#x00BB; should be &#x00AB;&#x00A0;Nik is the number of times
%(...)&#x00A0;&#x00BB;.
\end{document}


% LocalWords:  sidestepped sampleable Lennard alanine dipeptide Metropolized FL
% LocalWords:  Indep nonnegligible Bussi Roitberg Gervasio Laio Parrinello Chem
% LocalWords:  metadynamics Soc CVs al Gaussians CV PRL AE Simmerling Lett AB
% LocalWords:  BB ref Nij Nik alchemical
