{'fraction_non_alphanumeric': 0.03502688537164132, 'fraction_numerical': 0.013135082014365498, 'mean_word_length': 4.350403966998675, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 40, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computing reaction rates in bio-molecular systems using discrete macro-states

Computing reaction rates in biomolecular systems is a common goal of molecular dynamics simulations. The reactions considered often involve conformational changes in the molecule, either changes in the structure of a protein or the relative position of two molecules, for example when modeling the binding of a protein and ligand. Here we will consider the general problem of computing the rate of transfer from a subset A of the conformational space Omega to a subset B of Omega. It is assumed that A and B are associated with minimum energy basins and are long-lived states. Rates can be obtained using many different methods. We review some of the most popular approaches. We organize the different approaches roughly in chronological order and under four main categories: reactive flux, transition path sampling, conformation dynamics. The fourth class of methods, to which we do not give any specific name, in some sense attempts to combine features from transition path sampling and conformation dynamics. They include non-equilibrium umbrella sampling (Warmflash et al. [2007], Dickson et al. [2009b]), and weighted ensemble dynamics (Huber and Kim [1996]).


Introduction
Computing reaction rates in biomolecular systems is a common goal of molecular dynamics simulations. The reactions considered often involve conformational changes in the molecule, either changes in the structure of a protein or the relative position of two molecules, for example when modeling the binding of a protein and ligand. Here we will consider the general problem of computing the rate of transfer from a subset A of the conformational space Ω to a subset B ⊂ Ω. It is assumed that A and B are associated with minimum energy basins and are long-lived states. See Table 1 for the notations used in this paper.
Rates can be obtained using many different methods. In this paper we will review some of the most popular approaches. We organize the different approaches roughly in chronological order and under four main categories: reactive flux, transition path sampling, conformation dynamics. The fourth class of methods, to which we do not give any specific name, in some sense attempts to combine features from transition path sampling and conformation dynamics. They include non-equilibrium umbrella sampling (Warmflash et al. [2007], Dickson et al. [2009b]), and weighted ensemble dynamics (Huber and Kim [1996]).
Reactive flux. We start with methods that were derived around 1930s (Marcelin [1915], Eyring and Polanyi [1931], Eyring [1935], Horiuti [1938], Wigner [1938]), were revisited later by, for example, Chandler [1978], and are based on the concept of reactive flux. In these methods, the rate is derived from the free energy, and it is assumed that it is controlled by the flux at a saddle point at the top of the energy barrier separating A and B. The advantages of this approach is that it involves quantities which are relatively easy to calculate. However it makes relatively strong assumptions about the system, and in practice assumes that a lot of information is already available regarding the transition mechanism and important pathways between A and B. We relate some of the results to Kramers' method, which applies to systems modeled using Langevin dynamics and overdamped dynamics (Gardiner [1997], Hänggi and Borkovec [1990]).
Transition path sampling. Many of the ideas developed in that context were used to develop another class of methods based on sampling transition pathways between A and B , Dellago et al. [1998). From the ensemble of pathways, rates and other properties can be obtained. The advantage of some of these approaches is that they do not require determining the saddle point separating A and B, and they apply to more general situations, for example when multiple pathways contribute to the rate. See Vanden-Eijnden [2010] for a discussion of transition-path theory, which proposes a mathematical framework to study transition pathways and the associated probability current. To address shortcomings of some of these approaches, other methods were pursued along similar lines, including transition interface sampling (van Erp et al. [2003]), and forward flux sampling (Allen et al. [2005(Allen et al. [ , 2006a). We will present in the same category the milestoning technique, Faradjian and Elber [2004], which although different in spirit, shares some similarities with transition interface sampling. This approach requires that the system "loses" memory when moving from a milestone to the next, for example by ensuring that the milestones are sufficiently separated from one another. Recent advances will be presented based on the work of Vanden-Eijnden et al. [2008], who introduced the concept of optimal milestones using the committor function. Using these milestones, the rate can be obtained exactly (in the absence of statistical errors), even when the milestones are close to one another. Májek and Elber [2010] attempted to define the milestones in a way that would be computationally more general and advantageous compared to the original formulation.
Conformation dynamics. This is a large class of methods that can be traced back to Deuflhard and Schütte (Deuflhard et al. [1996], Schütte and Huisinga [2003]), and are based on the concept of metastable states and transfer operator (or transition matrix). Broadly speaking, Ω is decomposed into metastable sets, which are sets that are long-lived and in which the system gets trapped. Then a transition matrix P ij (τ ) is defined using the probability to reach a metastable set j if one starts a trajectory of length τ (the lag-time) in set i. The analysis of the eigenvalues lead to the concept of Perron cluster. From the eigenvectors and eigenvalues, one can derive the rate and other kinetic information.
Although derived apparently independently and at a later date, some groups started exploring how one could model molecular systems using Markov state models, a well-known theory but which has been only (relatively) recently applied to modeling bio-molecular systems. See Singhal et al. [2004], Swope et al. [2004a,b], Chodera et al. [2007Chodera et al. [ , 2006, Noé et al. [2007]. See related work by Shalloway [1996]. Many of the theory and results for Markov state models can be found in the literature on conformation dynamics. In some sense, Markov state models can be viewed as a practical implementation of conformation dynamics, that attacks the high-dimensionality of Ω by subdividing the space into a "small" number of cells, also called macro-states.
Specifically, the transition matrix P ij (τ ) giving the probability of reaching macro-state j when starting from state i after some lag time τ is used to compute the rate. Its first eigenvalue is equal to 1 and corresponds to the equilibrium distribution. The second eigenvalue is very close to 1 and can be used to estimate the relaxation rate, as well as the forward (reactant to product) and backward (product to reactant) rates. An important issue is the effect of the lag time τ . At short times, non-Markovian effects, or memory, are present, that is the Markov state model is not accurate and the estimated rate suffers from a systematic bias.
Reactive trajectory sampling. The last class of methods groups two separate approaches that in some sense combine ideas from transition path sampling and a subdivision of space similar to Markov state models. One such method, called weighted ensemble Brownian dynamics, originates in Huber and Kim [1996]. Although this paper is similar in spirit to transition interface sampling or milestoning, it can be easily extended to a general partitioning of space, using for example Voronoi cells. This is an important extension since, as a result, the method remains efficient in cases where multiple pathways contribute to the rate or when the most important pathway is not known. This approach leads to a sampling of transition pathways between A and B and therefore does not rely on the Markovian assumption made in Markov state models. However like Markov state models, the efficiency of the sampling is improved by partitioning space into macro-states. A large number of walkers (simulations) are run in each macro-state. In order to maintain the population of walkers in each macro-state, a procedure was created to kill walkers in macro-states that are too crowded, and to split walkers when the number of walkers becomes too low. This method was recently revisited by ,  who showed how the original approach could be extended.
The technique of non-equilibrium umbrella sampling of Warmflash et al. [2007], Dickson et al. [2009b] is similar in spirit. It applies, like weighted ensemble Brownian dynamics, to non-equilibrium systems and systems with memory (the Markovian approximation is not required) and uses a partitioning of Ω into macro-states. In each macro-state, a large number of walkers are simulated. Each time a walker attempts to leave a macro-state, its position (and velocity if needed) is recorded. Then in order to restart (continue) the simulation, a random position is chosen from the set of walkers who attempted to enter this macro-state from other macro-states.
In the second half of the paper we will discuss in more details weighted ensemble Brownian dynamics, renamed Reactive Trajectory Sampling (RTS) to reflect its broader application. The original method, Huber and Kim [1996], involves a procedure to split and kill walkers. We will revisit this method and propose an optimal procedure which leads to walkers with identical weights in each macro-state, a strategy which minimizes the statistical errors. We will discuss how the choice of macro-states affects the statistical errors and what the optimal choice is. We will present a new "multi-colored" algorithm that allows computing multiple rates (eigenvalues of the transfer operator from conformation dynamics) and accelerate convergence. This approach shares some similarities with the technique of core sets, .
We will propose a novel error analysis of Markov state models, by considering the sensitivity of the eigenvalues to perturbations in the transition matrix. This will lead to estimates of the systematic errors (non-Markovity) and statistical errors, and their dependence on the lag-time τ (length of trajectories used to calculate the transition matrix). The choice of macro-states influences the decay of the non-Markovian effects; we will discuss what the optimal choice is. This optimal choice is, as can be expected, difficult to realize in practice, but this provides nonetheless a guideline to improve and refine an initial guess. We will make an argument showing that in the general case, statistical errors increase with τ , showing that an optimal tradeoff must be found between memory effects (small τ ) and statistical errors (large τ ). The reader is referred to Sarich et al. [2010], Prinz et al. [2011], Hinrichs and Pande [2007], Singhal and Pande [2005], Metzner et al. [2009a] for papers that discuss the numerical errors in Markov state models. Some numerical results on simple benchmark problems in 1D and 2D are given at the end to illustrate the numerical properties of Markov state models and RTS.
Since there is a large number of methods to choose from, with different strengths and weaknesses, we attempt to summarize their main features for the purpose of comparing these methods together. We considered three axes in our evaluation: generality, computational cost, and parallel scalability. Each characteristic is ranked low, medium, or high. Generality relates to the number of approximations or assumptions that are required by the method to be accurate. For example, the method of reactive flux focusses its analysis on the transition region, typically a saddle point. When the assumptions are satisfied the calculation may be very accurate. However, in some instances the prediction may not be satisfactory and need to be improved for example using a method based on transition path sampling. The computational cost should also be taken as a general guideline since it will vary tremendously depending on the system. However some overall conclusions can be made regarding computational cost. Typically the situation is that generality is traded for computational cost. Finally since all large scale calculations require a parallel computer (multicore processors, graphics processing units, and parallel clusters), we also rank methods according to the amount of parallelism they offer, although in all cases, the amount of parallelism is very large and scalability is typically not an issue. Codes that run on slow networks (grid computing, cloud computing) will be more sensitive to these issues. Finally, the information in the list below should be taken merely as a guideline since all conclusions are in general strongly system and implementation dependent.

Reactive flux
Generality: low. These approaches are typically the ones that require the greatest amount of knowledge about the system and relatively strict assumptions, in particular regarding the energy landscape near the transition region. Computational cost: low. A free energy calculation is required. Although this can be difficult in some instances, this type of calculation is typically easier than with the other methods. Parallel scalability: medium. Many methods are available to calculate free energy and are quite scalable. In most cases, it may become difficult to increase the number of processors compared to the other techniques described here.
Transition path sampling Generality: high. These approaches are among the most "direct" and require little or no assumption. Computational cost: high. As a result of being very general, they lead to extensive sampling and typically long simulation times. Parallel scalability: high. In most cases, one can sample paths independently.
Transition interface sampling, forward flux sampling Generality: high. These approaches also require little or not assumption. Forward flux sampling in addition only requires the ability to run forward simulations making them applicable in situations where transition interface sampling fails. Computational cost: high. A large number of pathways must be sampled. Forward flux sampling may in some cases converge a little slower, in particular when the sampling of the initial interfaces is poor or turns out to be insufficient for later interfaces. These methods typically perform better than transition path sampling. Parallel scalability: medium. Processing the interfaces is sequential, making the method less parallel.
Milestoning Generality: medium. Strong assumptions are required regarding loss of memory between milestones. As discussed in the main text, the original method has been extended and made more general. Computational cost: medium. Only short pathways between adjacent milestones are required making this approach less expensive than the two previous methods. Parallel scalability: high. The milestones can be processed independently.
Markov state models Generality: medium. The Markovian assumption must apply, which depends on the choice of Markov (or macro) states and the lag time between observations. Computational cost: medium. The sampling is mostly local within each macro state, leading to efficient sampling. Parallel scalability: high. Macro states are processed independently and require running a large number of short trajectories (in some implementations).
Weighted ensemble Brownian dynamics and non-equilibrium umbrella sampling. Generality: high. The Markovian assumption is not required. The level of generality is similar to transition path sampling and related methods. The convergence of weighted ensemble Brownian dynamics is relatively easy to monitor. Computational cost: high. The fact that the Markovian assumption does not apply typically leads to higher computational cost compared to Markov state models. Parallel scalability: high. It is similar to Markov state models, although some communication is required to update the weight of walkers. This involves a small amount of communication, but one that occurs at regular intervals.
Weighted ensemble Brownian dynamics is in many respects similar to transition path sampling and differs "mostly" in the technique used to generate paths joining A and B and enhance the sampling. In this paper, the method of weighted ensemble Brownian dynamics will be referred to as reactive trajectory sampling to indicate that more general formulations have been created since the original paper of Huber and Kim [1996].
Many theoretical results for these methods have been proved in the context of Langevin dynamics or Brownian dynamics. Results are often derived in the context of Brownian dynamics (over-damped dynamics). Extensions to Langevin are in most cases possible, although the proofs become more technical. In some cases, a method may depend only on some stochastic process with minimal assumptions. The theoretical derivations at the end of this paper are done in the context of Brownian dynamics but results can be extended to Langevin dynamics. Extensions to Newtonian dynamics (deterministic) are much more difficult and in most cases these extensions do not exist yet.
Lag-time in Markov state models P ij (τ ) probability to be in state j when starting from i after time τ ξ reaction coordinate or order parameter that monotonically increase from A to B (ξ 1 , . . . , ξ p ) set of generalized coordinates A(ξ) free energy D diffusion tensor χ A characteristic function of some set, A C(t) conditional probability to find the system in B at t provided it was in A at time 0 x(T ) a discrete trajectory in Ω P(x(T )) probability density function of trajectories P AB [x(T )] probability density function for the transition path ensemble S i hypersurface in Ω; used in transition interface sampling, forward flux sampling, etc. V i cell or macro-state in Markov state models and related methods In many cases S i is the set of points such that ξ(x) = ξ i Φ A,1 flux from region A to S 1 P A (ξ i+1 |ξ i ) for trajectories coming from A, probability to reach S i+1 starting from S i π(x) committor function µ i eigenvalue of the transition matrix λ i λ i = −(ln µ i )/τ ; they are often an approximation of the eigenvalues of the Fokker-Planck equation P( ) used to denote the probability of some event happening X n , X t homogeneous Markov process indexed by n or t N ij number of observed crossings from macro-state i to j T i length of simulation in macro-state i W i statistical weight of macro-state i w i statistical weight of walkers in weighted ensemble Brownian dynamics or reactive trajectory sampling ρ(x, t|x 0 , 0) probability to be at x at time t if the system was at x 0 at time 0 ρ k eigenfunction of the forward Fokker-Planck equation ψ k eigenfunction of the backward Fokker-Planck equation λ k corresponding eigenvalues 2 Transition path sampling

Reactive flux and transition state theory
We start the discussion with the method of reactive flux which is a long standing approach to computing reaction rates. The idea goes back to Marcelin [1915], Eyring and Polanyi [1931], Eyring [1935], Horiuti [1938], Wigner [1938], who developed the initial theory of chemical reaction kinetics. We outline the main ingredients in this type of approach. We assume that region A is a subset of the conformational space of the molecular system and that it represents in the system in its reactant state. Similarly B denotes the region defining the product states. Analytical approximation for the rate can be obtained if one assumes that a coordinate ξ can be defined which describes the reaction, a reaction coordinate. It is assumed that when ξ = 0 the system is in A and when ξ = 1 the system is in B. The value ξ = ξ * corresponds to the transition region or barrier between A and B.
We define the characteristic function χ A (resp. B) which is 1 in the set A and 0 outside. Then using these functions, we can express the conditional probability to find the system in state B at time t provided it was in A at time 0: Brackets are used to denote a statistical average. Regions A and B are separated by a transition region and the rate is determined by the rate at which this transition or barrier is crossed. At the molecular scale, there is some correlation time τ mol associated with this crossing. That is for times larger than τ mol , the system has forgotten how it went from A to B. Then for times t between τ mol and the reaction time τ rxn , τ mol < t τ rxn , the time derivative of C(t), called the reactive flux, reaches a plateau (Chandler [1978]), The symbol˙denotes a time derivative. Using transition state theory (TST), under the assumption that the recrossing of the barrier between A and B can be neglected, one can derive an expression for k AB using Eq. (1) (Chandler [1978(Chandler [ , 1987): where A(ξ) is the free energy, and ξ=ξ * denotes an ensemble average with ξ constrained at ξ * . This approach has some drawbacks. It always overestimates the rate. It requires a good reaction coordinate and a precise determination of the free energy maximum to locate the barrier. Nevertheless the method is computationally efficient and involves only quantities that can be computed with relatively low computational cost. Among the many methods to calculate the free energy (in this context the potential of mean force), see for example Lelièvre et al. [2010], Chipot and Pohorille [2007], Darve and Pohorille [2001], Henin and Chipot [2004], Rodríguez-Gómez et al. [2004], Lelièvre et al. [2007], Darve et al. [2008]. Related approaches include Kramers' rate theory (Gardiner [1997], Hänggi and Borkovec [1990]), which was developed in the context of Langevin equations and overdamped dynamics. There are many connections between transition state theory and Kramers' theory, Hänggi and Borkovec [1990]. In particular Kramers' rate can be related to the "simple" TST rate through: In this expression the potential at the transition point is assumed to be locally quadratic with stiffness ω 2 bar = −m −1 U (x bar ) (m is the mass of the particle in a 1D model), and λ + is a function of the friction in the Langevin model and ω bar . It can be shown that k Kramers is equal to the multidimensional TST rate for a heat bath describing strict Ohmic friction, Hänggi and Borkovec [1990], pp. 268 & 272. As the friction in the Langevin model goes to zero λ + → ω bar and k Kramers → k TST . Moreover we always have k Kramers < k TST . The rate k Kramers is itself an upper bound on the true rate given by where θ is the Heaviside function. This basic approach using TST has been improved in many ways including the use of harmonic approximations to model the minimum energy basins and transition region; see Dellago and Bolhuis [2009], Gardiner [1997], Hänggi and Borkovec [1990]. In variational TST, one attempts to improve the predicted rate by finding a dividing surface between A and B that minimizes the rate, see for example Truhlar [1984], Tucker [1995]. The overestimation of the rate by TST is a result of neglecting the re-crossing of the dividing surface. Several authors have proposed corrections to the basic TST approach to account for these effects. See for example Bennett [1977], Chandler [1978].

Transition path sampling
The method of reactive flux is attractive as its computational cost is often tractable. The primary calculation is obtaining the free energy profile along the reaction coordinate ξ. The accuracy is very dependent on the choice of coordinate. Specifically, the transmission coefficient, which measures the amount of re-crossing taking place, is dependent on the definition of ξ and a low transmission coefficient leads to inaccuracies or inefficiencies. Transition path sampling were in part proposed to alleviate the need to define this coordinate as they rely primarily on sampling trajectories going from A to B with no knowledge of ξ required. However we will see that later versions, again, are dependent on some knowledge of ξ.
Transition path sampling , Dellago et al. [1998) is a Monte-Carlo method that allows sampling the ensemble of trajectories. For each discrete trajectory one can define a probability P(x(T )) to observe such a trajectory. Then the transition path ensemble defines a probability density in the space of trajectories that is non-zero only for trajectories that connect A and B, and therefore its probability density function is defined as: where Z AB is the appropriate normalization factor. Trajectories in this ensemble can be generated using a Monte-Carlo procedure. A popular method is the so-called shooting method, in which a point is selected along the trajectory, the momentum is modified and a new trajectory is generated using a forward and backward time integration. If a stochastic dynamics is used then no perturbation is necessary since the random number generator will lead to a different trajectory. Then the appropriate acceptance probability is used (Metropolis-Hastings algorithm) to accept or reject this new trajectory. See Fig.1. Trajectories that do not start in A and end up in B are always rejected. This procedure can be improved by adding extra moves like shifting moves or path reversal moves, .

A B
Point at which perturbation in momentum is applied Figure 1: Shooting method of transition path sampling. The red curve is the old trajectory. The blue curve is obtained by perturbing the momentum at some point and regenerating a trajectory by using a forward and backward time integrator.
Many thermodynamic and kinetic properties can be determined from transition path sampling. As before the rate is related to C(t) through k AB (t) =Ċ(t) and (Dellago and Bolhuis [2009], ) where Dx(t) is an integration over all possible paths (precisely, all possible points along the approximation at discrete time steps of the continuous trajectory).
There is a numerical difficulty in computing since for most trajectories of length t, χ A (x 0 ) = 0 or χ B (x t ) = 0. A common procedure to make this tractable is to use umbrella sampling and an order parameter ξ which approximates the reaction coordinate. Denote P A (ξ , t) the probability that a trajectory started in A is such that ξ = ξ at time t, then (Dellago and Bolhuis [2009]): where ξ B min ≤ ξ ≤ ξ B max defines region B. The quantity P A (ξ, t) will be difficult to obtain for values of ξ that are close to 1/2 (barrier) or beyond. The method of umbrella sampling (Torrie and Valleau [1977]) can then be conveniently used. For a window W i in the interval [0, 1] define: Then each P Wi A (ξ, t) can be efficiently computed for small enough windows, while P A (ξ, t) over the entire [0, 1] interval is obtained by patching together the different P Wi A (ξ, t) and matching the curves to produce a single profile.
Transition path theory. Transition path theory (TPT) is a theoretical framework to study transition path ensembles. It considers so-called reactive trajectories, which are the trajectories sampled by transitionpath sampling. The theory is derived in the context of Langevin and over-damped dynamics. TPT provides definitions for the probability density of reactive trajectories. Most results involve the forward q + or backward q − committor functions. The forward committor function for example is defined as the probability that, starting from some point outside of A ∪ B, the system first reaches B before A, Vanden-Eijnden [2010]. As an example, the probability density of reactive trajectories is defined as: where ρ(x) is the equilibrium distribution for the process. More importantly expressions are provided for the probability current of reactive trajectories. This allows identifying important transition pathways and transition tubes. The example of the maze is helpful. The path in the maze with the largest reactive flux corresponds to the shortest path from A to B. Actual transition pathways take many detours and visit many dead ends before back-tracking and ending up in B. Analyzing the reactive flux allows bypassing these detours and take a straight route to B. In practice this allows identifying important reaction mechanisms joining A and B. Expressions are also given for the reaction rate, basically by integrating the probability flux over a subdividing surface. Metzner et al. [2009b] presents an application to discrete Markov processes. In this case the discrete probability current for reactive trajectories is given by: where l ij is the infinitesimal generator (

Transition interface sampling
One issue in the previous approach is the fact that the trajectories need to have a fixed length t. See  (Section IV D "A Convenient Factorization") for an algorithm that relaxes this requirement. This length must be chosen carefully. A short time will lead to inaccuracies while a long time leads to a larger computational cost to update the trajectories. The method of transition interface sampling (TIS) introduces a novel concept, somewhat related to the previous algorithm with umbrella sampling, in which n + 1 multi-dimensional surface in the space outside of A and B are defined. These surface are such that they can be used roughly to measure the progress of the reaction. For example, for some order parameter ξ and choosing a value ξ i , we can define a surface S i by the equation ξ(x) = ξ i . We assume that ξ 0 < ξ 1 < . . . < ξ n ; S 0 is the boundary of A and S n is the boundary of B. Several methods (transition interface sampling and forward flux sampling) start from a reinterpretation of the rate as: The first term, Φ A,1 , is the flux of trajectories going out of A and crossing the first surface S 1 . This quantity can be obtained by running a long trajectory (ignoring parts of the trajectory that last visited B) and counting the number of times the surface S 1 is crossed, per unit time. Only positive crossing are counted (that is moving away from A towards B) and re-crossings are ignored until A is entered again (that is crossings for which the previous crossing was S 1 , not A, are ignored). The second term, P A (ξ n |ξ 1 ) is the probability to reach S n , assuming that the trajectory has crossed S 1 and that it does not cross A before S n . See van Erp et al. [2003]. The probability P A (ξ n |ξ 1 ) is difficult to calculate since it is very small. This can be remedied using the intermediate surfaces S i , 2 ≤ i ≤ n − 1. By construction, the function ξ being continuous and since ξ 0 < ξ 1 < . . . < ξ n , it is not possible to last come from A and cross ξ i without having crossed first ξ i−1 . Using this result, it is possible to show that (van Erp et al. [2003]): where P A (ξ i+1 |ξ i ) is the probability to cross S i+1 before A assuming the trajectory had previously crossed A then S i . Although P A (ξ n |ξ 1 ) can be very small, the quantities P A (ξ i+1 |ξ i ) are much larger and can be reliably obtained through direct sampling. The procedure to calculate P A (ξ i+1 |ξ i ) is similar to the procedure for TPS. It is illustrated in Fig. 2a.

Forward flux sampling
The forward flux sampling method (FFS) was conceived by Allen et al. [2005Allen et al. [ , 2006a]. An earlier paper by Harvey and Gabb [1993] bears some conceptual similarities. In TIS, one generates new trajectories by perturbing a point and then integrating forward and backward. In the forward flux sampling method, only forward integration is used. This can be essential in cases where backward integration is not possible. For example in non-equilibrium systems, the lack of detailed balance and absence of time-reversal symmetry means that TIS (or milestoning which will be described in the next section) is not applicable. FFS is one of the few methods applicable to such systems. FFS uses the same basic framework as TIS expressing the rate as FFS start by calculating P A (ξ 2 |ξ 1 ) using trajectories initiated from S 1 . Some of these trajectories may fail to reach S 2 (i.e., reach A before S 2 ) while others will reach S 2 successfully (before reaching A). For those that reach S 2 , the first hitting point (first point where the trajectory crosses S 2 ) is saved. Those points are then used to calculate the next conditional probability, P A (ξ 3 |ξ 2 ). FFS therefore only requires forward integration of trajectories. The starting points on S i are in some sense fixed and are produced solely as a result of the sampling during the calculation of P A (ξ i |ξ i−1 ). Typically the computation for P A (ξ i |ξ i−1 ) is continued until satisfactory accuracy is achieved and enough points have been generated on S i . We note that the forward flux approach not only yields the rate constant, but also the complete transition pathways, which can be reconstructed by "gluing" the successful trajectory segments together. The method is depicted in Fig. 2b. A limitation of this method is that the accuracy at later interfaces depends on the sampling at earlier interfaces. For example if the first interface is relatively poorly sampled the error will propagate throughout the next interfaces. In addition it is very possible that initial trajectories with low probabilities, in the end, make large contribution to the flux. Such a situation would lead to a large standard deviation and statistical errors.
Example applications are discussed for example in Valeriani et al. [2007], Borrero and Escobedo [2007], Allen et al. [2005]. Allen et al. [2006b] describe FFS along with two other methods, the branched growth method and the Rosenbluth method. Allen et al. [2006a] proposes an analysis of the efficiency of these methods.

Milestoning
Similar to the previous method, the milestoning technique of Faradjian and Elber [2004], West et al. [2007] is based on a set of separating hyper-surfaces that are used to measure the progress of the system from A to B. The advantage of this method is that trajectories need only to be run from one interface S i (or milestone) to the next S i+1 or previous interface S i−1 . In contrast, the previous methods required running trajectories until they reach the next milestone or set A. Milestoning is inherently more parallel since each interface S i can be processed in parallel with the others, whereas the other methods require a sequential processing of the interfaces. A drawback is that the interfaces need to be sufficiently separated so that we can assume that the system loses memory in the time it takes to reach the next interface.
The milestoning method proceeds by initiating trajectories at S i , using the equilibrium distribution. Then it records the time it takes to reach S i−1 or S i+1 . This leads to two time distribution functions Fig. 2c. Then we can calculate two functions: P s (t), which is the probability of being at milestone s (that is the last interface that was crossed was S s ), and Q s (t), the probability to transition to milestone s at time t (cross S s at t). Then these two functions satisfy the following set of integral equations: In these equations η s is the initial milestone probability distribution. The first equation is saying that in order to be at s the system must first transition to s [Q s (t )] and then stay at s until time t > t . The second equation says that to reach s one must first get to s − 1 and then transition from s − 1 to s (and similarly with s + 1). The first term η s δ(t − 0 + ) accounts for the initial distribution at t = 0 (δ is the Dirac delta function). From these equations, the free energy and reaction rate can be obtained. This approach is accurate once we assume that the system loses memory between milestones. With this, it becomes justified to independently generate initial conditions on each milestone following the Boltzmann distribution. These assumptions allow running independent calculations at each milestone and make it possible to run trajectories that stop as soon as the previous or next milestone is reached.
See a discussion of this approach in Elber [2005] as well as an application example to an allosteric transition with deoxy Scapharca hemoglobin in Elber [2007].

Milestoning using optimal milestones
The method of milestoning was recently revisited, see e.g., Vanden-Eijnden et al. [2008]. In this paper, it is shown that the assumption that the milestones need to be sufficiently far apart is not required provided that the interfaces S i are iso-surfaces of the committor function. The committor function, denoted π(x), is the probability to reach B before A starting from x. In Vanden-Eijnden and Venturoli [2009a], the milestoning (a) Transition interface sampling algorithm. Calculation of P A (ξ 3 |ξ 2 ). The sampling is similar to TPS. The main difference is that the only paths that are considered are those that last come from A, cross S 2 and then cross A or S 3 . A new path is generated by perturbing a given path and then integrating forward or backward.
The integration can be stopped as soon as A or S 3 are reached. The new path is then accepted only if there is a segment crossing S 2 that last come from A. Green trajectory: it has not reached the surface S 2 and therefore is discarded. Red trajectory: it has reached S 3 before A and therefore counts as 1. Blue trajectory: it has reached A before S 3 and therefore counts as 0. (c) Milestoning algorithm. Trajectories are initiated from S 2 using the equilibrium Boltzmann distributions. The exit times are recorded when the system reaches S 3 (red trajectories, function K + 2 (t)) or S 1 (blue trajectories, function K − 2 (t)). method is extended to the case of a tessellation of the conformational space using Voronoi cells V i . See Fig. 3. See a more complete discussion about Voronoi cells in Section 3.3. In a more recent paper, Májek and Elber [2010] modify the original milestoning method using an approach that does not require a reaction coordinate (in that sense similar to Vanden-Eijnden and Venturoli [2009a]). However, instead of considering the committor function to guarantee the accuracy of the method, this approach focusses in guaranteeing a minimum separation between the milestones. This leads to greater accuracy since the assumption that memory is lost between crossing milestones is automatically satisfied. We now review some of the mathematical underpinnings discussed in Vanden-Eijnden et al. [2008]. In particular we recall the main arguments to establish that the use of the committor function π(x) to define the milestones leads to an exact rate prediction with milestoning, irrespective of whether the milestones are close or well separated. For this reason these milestones are called optimal milestones.
In Vanden-Eijnden et al. [2008], one is concerned with computing the mean passage time τ AB from A to B.
Consider for example the case of over-damped dynamics and the associated Fokker-Planck equation. The eigenvalues of the Fokker-Planck equations are denoted λ i . We assume the λ i are ranked by magnitude so that λ 1 = 0 < λ 2 < λ 3 < . . .. If there is a single eigenvalue λ 2 close to 0 (λ 2 λ 3 ), then the mean passage time is related to λ . Milestones, denoted by S i , are defined as hypersurfaces in Ω. The mean time to go from A to B can be obtained from the mean times to go from a milestone to another milestone. This has important consequences in terms of which assumptions need to be made to get an exact rate (in the absence of statistical errors).
The key property is the following one. Let us assume we consider all the trajectories that go from milestone S i to S j (j = i + 1 or j = i − 1). Specifically we initialize trajectories on milestone S i with density: where dσ i (x) is the surface element on milestone S i . Then we can follow all the trajectories and record where they hit milestone S j (we discard trajectories that hit other milestones first). The density of points on S j is called the density of first hitting points.
In Vanden-Eijnden et al. [2008], it was shown that, if the milestones {S k } are iso-surfaces of π (optimal milestones), then the density on S j is given by: In particular this density is independent of the fact that the trajectories were started from S i . This can be proved using the definition of the committor function and the forward Chapman-Kolmogorov equation.
The implication is that if we want to calculate the mean time to go from i to j we do not need to consider where the system is coming from. It is sufficient to initiate trajectories on S i with the density given above and calculate the average time required to reach S j (assuming this is the next milestone crossed).
We now discuss the calculation of the mean passage times and show that the property above is essential to derive expressions for the exact rate. Let us consider an absorbing boundary condition on the milestone, denoted S n , that surrounds region B (cemetery milestone Vanden-Eijnden et al. [2008]). We assume that trajectories are initialized on S i with a probability density given by Eq. (17). We define T i the mean length of a trajectory that start from S i and goes to S n (mean exit time) . We denote: p ij the probability that S j is crossed after S i , and τ i the mean time to hit any other milestone S j , j = i. Using the result above regarding the first hitting point density, we have: This is true because the mean escape time for some S j , j = i, is the same irrespective of whether the system comes from S i or not. This is a weaker property than saying that the times to go from S i to S j , t ij , and S j to S k , t jk , are independent. One can construct examples where t ij and t jk are strongly correlated, while the density of first hitting points at a milestone is independent of the last milestone that was crossed. This equation can be written in matrix form as: where p n ii = 0 and p n does not have the row or column corresponding to S n . This equation requires computing p ij and τ i for all milestones. This can be further simplified by observing that p ij /τ i has a simple interpretation. A formal proof can be given but here we simply outline the main points. We consider a very long trajectory. We will not discuss this further but limits must be taken as the trajectory length goes to infinity. Then: Number of times the system crossed S j after S i Number of times the system crossed any milestone after last crossing S i (21) τ i = Total length of time during which the last milestone crossed was S i Number of times the system crossed any milestone after last crossing S i So that: p ij τ i = Number of times the system crossed S j after S i Total length of time during which the last milestone crossed was S i Assume that we now use a discrete integrator to integrate the dynamics (which may be Langevin or overdamped dynamics), with time step ∆t. We define: P ij = P(system crosses S j during the next step assuming it last crossed S i ) (24) P ii = P(system crosses no milestones during the next step assuming it last crossed S i ) Then from Eq. (23): if ∆t is small compared to the time required to go from a milestone to the next. This shows that: 1 The quantity ∆t −1 T i is the mean number of steps required to go from S i to S n . The advantage of this equation is that it requires computing P ij only, which can be relatively easily computed. Assume that we have defined a partition of the conformational space Ω into cells V i such that the milestones S i form the boundary of these cells. See Fig. 3. Then the approach requires simply running independent simulations in all the cells V i . For this, we need to use boundary conditions such that the system remains in the cell it started in, during the simulation. Let us assume that the trajectory hits a cell boundary with velocityẋ. Then we know from the equilibrium probability density that there is another trajectory in the past or future, with a one-to-one mapping, which re-enters the cell through the same point. Its velocity can be chosen equal to −ẋ (using the fact that the equilibrium probability density is even with respect to the momenta), or obtained using a hard wall reflection (now using the fact that the reflection conserves the equilibrium probability density) witḣ x re-entering particle =ẋ − 2 (ẋ · n) n and n is the normal to the interface.
With this approach, one can generate a large number of samples in each cell, from which we can estimate P ij : where N ij is the number of times the system was found to cross S j after S i , a is the cell bordered by S i and S j , while b is the cell on the other side of S i , N a i (resp. N b i ) is the number of steps for which the last milestone crossed was S i in cell V a (resp. V b ), and n a , n b are the number of steps computed in each cell. This is basically a direct calculation of Eq. (23). Note that by construction, the density of first hitting points on the milestone is the exact one so that, up to statistical errors, Eq. (29) is exact.
The advantage of this approach is its efficiency and the fact that the accuracy is more or less independent of the energy barrier between A and B. It does not have any systematic error unlike the previous approaches. The main drawback is the requirement that the optimal milestones are iso-surfaces of the committor function, which again is difficult to realize in practice.
This requirement was relaxed in Májek and Elber [2010].
3 Conformation dynamics and Markov state models

Conformation dynamics
All the methods discussed above attempt in a sense to do a direct calculation of the rate, either by computing a reactive flux at the transition barrier or by sampling reactive trajectories going from A to B. We now discuss another class of methods that also attempt to calculate the free energy and the rate but, indirectly, by calculating the rate of transition between metastable basins. If one derives a statistical model of the system in terms of hops or transition between states, then an eigenvalue analysis can be used to calculate reaction rates, metastable states, and extract many relevant kinetic and thermodynamic information.
The idea goes back to Deuflhard and Schütte, who realized that computing time averages of physical observables or relaxation times of conformational changes (using molecular dynamics for example) was largely determined by the existence and properties of invariant sets, called metastable sets. These are by definition sets (subsets of Ω, the conformational space of the molecular system) such that the system stays trapped in these sets for extended periods of times and with very rare transitions between sets. This has led to the conformation dynamics approach which aims at identifying these sets, and computing the transition rates between these sets. The first paper goes back to Deuflhard et al. [1996], although the term "essential dynamics" can be found in Amadei and Linssen [1993]; Grubmüller and Tavan [1994] had introduced the concept of conformational substates in 1994. Deuflhard [2003], Schütte and Huisinga [2003] provided some surveys on this topic. Central to this model is the concept of transfer operator and the study of its eigenvectors and eigenvalues. Relevant reaction rates (as well as mean passage times, mean exit times, . . . ) can then be extracted from these eigenvalues. We provide a brief account of the key mathematical objects. The definitions can be made in a relatively general context (Schütte and Huisinga [2003]). We consider a stochastic transition kernel p(t, x, A) such that: where the family {X t } is a homogeneous Markov process indexed by a time variable t; p(t, x, A) is therefore the probability that a Markov process started at x is in A after a time span t. This allows defining the Perron-Frobenius operator P t (propagator or forward transfer operator). See Schütte and Huisinga [2003], Schütte et al. [2001] for a definition that does not assume that the transition kernel p(t, x, y) is continuous with respect to a probability measure µ(dx). Then: The density ρ is assumed to be invariant. In a similar fashion, ensemble transition probabilities can be defined as:

Perron cluster cluster analysis
The metastable sets alluded to earlier are defined as sets that are almost invariant under the Markov process. With our definition, the set C is almost invariant if: These almost invariant sets can also be identified by considering the eigenvalues and eigenvectors of P t . For example, the density ρ, being invariant, satisfies P t χ Ω = χ Ω ,where Ω is the space of all conformations of the molecule and χ Ω is the characteristic function of Ω (in that case simply the function 1). The eigenvalue 1 therefore corresponds to the equilibrium distribution ρ.
Eigenvalues close to one form the so-called Perron cluster. The largest of these eigenvalues (not equal to 1) can then be associated with the slowest rate in the system. See the thesis of Huisinga [2001] for an extensive discussion of transfer operators and metastability. Huisinga and Schmidt [2006] provides some upper and lower bounds on these eigenvalues. There is a significant body of literature on the identification of these sets, in particular using an approach called Perron cluster cluster analysis (PCCA). The premise is that if one has uncoupled Markov chains (that is P t is a block diagonal matrix assuming discrete states) then, assuming k separate aggregates or clusters, we will have k eigenvectors with eigenvalue 1. Each eigenvector is constant over each cluster and changes sign (positive, negative or 0) between clusters ). In PCCA, the sign has been used to identify these clusters or aggregates.
In a real application however, we are not dealing with uncoupled Markov chains but rather weakly coupled chains, resulting in a single eigenvector with eigenvalue 1 and a cluster of eigenvalues near 1, Schütte and Huisinga [2000b]. The identification of the clusters based on the sign structure ) is then more difficult as the sign change is more progressive with a smooth change of value across the transition region, and the determination of the sign is more difficult as the eigenvectors may assume very small values whose sign (+1, −1, 0) is then difficult to determine.
One remedy to this is to recognize, as will be further discussed later on, that the sign of an eigenvector entry is not the right quantity to look at. For example the second eigenvector ρ 2 with eigenvalue 1 − ε is typically nearly constant with value say ρ A 2 > 0 in one cluster and is constant with value ρ B 2 < 0 in another cluster. The transition region is not associated with the point where ρ 2 is zero (or changes sign) but rather where ρ 2 assumes the value 1/2(ρ A 2 + ρ B 2 ) (midpoint between the two plateaus). This midpoint value can be used in a robust fashion to determine the cluster boundaries.
Another approach pioneered by Deuflhard (Deuflhard and Weber [2005], Weber and Rungsarityotin [2004]) called robust Perron cluster cluster analysis (PCCA+) proposes a fuzzy decomposition where instead of a strict partitioning into clusters one calculates of partition of unity iχ i (l) = 1 (at a discrete state indexed by l) where each functionχ i is called an almost characteristic function which smoothly transitions from 1 to 0 outside out a cluster. In effect, Weber and Rungsarityotin [2004] assign a grade of membership between 0 and 1 to each state (in a discrete setting). Therefore, each state l may correspond to different clusters with a different grade of membership, defined byχ i (l). This approach was used to study a relatively long polyalanine (Ala 8 and Ala 12 ) in Noé et al. [2007]. The backbone torsion rotamer pattern was used to define the microstates.
The thesis of Weber [2007] reviews these methods and discusses a meshless approach in which the membership functions are used to construct a Galerkin approximation of P t . The fact that the basis functions are non-orthogonal (in contrast with partition functions for example) leads to a non-diagonal mass matrix in the Galerkin formulation.
We mention an alternative method based on a singular value decomposition, see Fritzsche et al. [2007]. The main drawbacks of the approaches mentioned above (PCCA and variants) are the difficulty of identifying the Perron cluster if the transition matrix of the Markov chain has no significant spectral gaps; in addition, the calculation of the eigenvectors may be badly conditioned if the Perron cluster contains many eigenvalues very close to 1. The SVD approach attempts to mitigates some of these issues. Weber [2005, 2007]) have used this decomposition into metastable sets with PCCA+ to construct a coarse grained matrix that approximates the exact fine grained propagator. Limitations of the resulting propagator are discussed.
An application of this approach to systems described by a Langevin equation is given in Schütte and Huisinga [2000a]. An issue in this approach is the fact that the conformational space needs to be discretized appropriately in order to form a discretized approximation of P t (called transition matrix in Cordes et al. [2002]). This is difficult in practice since molecular systems live in high dimensional space. Cordes et al. [2002] propose to focus on dihedral angles to reduce the dimension and further improve their approach by considering a recursive decomposition in which space is first decomposed using the dihedral angle with the longest auto-correlation (which is shown to correlate with metastability). This leads to a first partitioning of space. Then, these metastable subspaces are further decomposed by applying the same strategy recursively (see Cordes et al. [2002]). PCCA is finally applied to the resulting coarse decomposition. Galliat et al. [2000] proposes an approach based on self-organized neural networks, also to attack this problem of dimensionality. Schultheis et al. [2005], Kloppenburg and Tavan [1997] approach the problem of dimensionality through the use of density-oriented discretizations that represent the probability density using a mixture of normal distributions.
The concept of metastability also allows viewing the Markov chains as a collection of rapidly mixing chains within certain parts of the state space, that are weakly coupled together. This gives rise to the concept of uncoupling wherein uncoupled Markov chains, which resemble the original chain, are formulated for each metastable states. Then these k chains can be coupled again by a k × k coupling matrix. The resulting system contains all the important information from the original chain. This is discussed in Meerbach et al. [2005], where in addition upper bounds are provided for the 2nd eigenvalue of the uncoupled Markov chains to establish that they are indeed rapidly mixing.
In a related approach to building coarse models without resorting to the transition matrix or transfer operator, Horenko et al. [2006] attempts a direct construction of a coarse grained model using a hidden Markov model with output given by stochastic differential equations (SDE). The model is then formulated as the combination of a set of SDEs (Langevin equations) and a rate matrix which determines how one jumps from an SDE to the next. See Horenko et al. [2006], Horenko and Schütte [2008]. The idea of hidden Markov model is reused in Fischer et al. [2007] but this time the output of a hidden Markov state is a probability density function in the observed variables (in this case torsion angles along the protein backbone). Von Mises output distributions are used since the observed variables are angles.
Although many clustering methods have been proposed, we mention the work of Coifman et al. [2008] who proposed to build a coarse grained representation based on the eigenvectors of the diffusion map. The algorithm is based on the definition of a weighted graph on the simulated points and the subsequent computation of the first few eigenvalues and eigenvectors of a random walk on this graph. Connections are made to the backward Fokker-Planck operator.
A related approach was developed by Shalloway and his group. See Orešič and Shalloway [1994], Church et al. [1996], Shalloway [1996], Ulitsky and Shalloway [1998], Korenblum and Shalloway [2003]. Church et al. [1999] review these methods. In Orešič and Shalloway [1994], Gaussian packets are used to characterize metastable basins and equations are provided to calculate and evolve packets. Packets are Gaussian functions in Ω. These packets can be numerically obtained by computing the average value of x locally in some metastable basin, and the variance of x: Orešič and Shalloway [1994] for the notation i (this is a suitable local averaging in basin i). Then the Gaussian packet has the form: In Shalloway [1996], the eigenfunction expansion of the Smoluchowski equation is used to optimize the macrostate expansion (using Gaussian packets) by satisfying a minimum uncertainty condition. Ulitsky and Shalloway [1998] builds macrostates based on a variational principle, from which the transition region (separating the macrostates) can be identified and analyzed. Connections with the variational transition state theory (Pollak in Fleming and Hänggi [1993], p. 5-41) are made.

Markov state models
Another broad class of methods are Markov state models (MSM). Markov chains have a long history, however their application to bio-molecular modeling and protein modeling is relatively recent and goes back to papers by Singhal et al. [2004], Swope et al. [2004a], and Swope et al. [2004b]. Two landmark papers are Chodera et al. [2007] and Chodera et al. [2006]. See  for a review and discussion of this model. This approach is closely related to conformation dynamics and many theoretical results from conformation dynamics directly apply to MSM. To build a discrete model of the continuous underlying conformational space Ω, the method starts by partitioning Ω into cells. The advantage is that this decomposition is often easier to obtain in practice than defining a reaction coordinate. In this respect this approach has some of the advantages of transition path sampling, which also does not require a precise knowledge of the reaction coordinate. This is a significant departure from TIS, FFS, and milestoning who rely primarily on a single order parameter that measures the progress of the reaction from A to B in an essentially sequential way (progress along a single dimensional variable or order parameter). MSM remains easy to construct even in the presence of multiple pathways with no obvious reaction coordinate. These macro-states can be constructed in different ways, for example from pathway sampling information (Singhal et al. [2004]) or by partitioning the Ramachandran map (Chodera et al. [2006]).  uses a decomposition based on torsion angles. In this paper, a method is also proposed to define generalized angle coordinates, in an attempt to reduce the number of coordinates in the problem. It is based on ideas from Amadei and Linssen [1993], and statistical analysis of circular data (Fisher [1993], Fisher and Lee [1983]). Chodera et al. [2007] proposes an automatic procedure to create macro-states using an iterative procedure and the k-medoids algorithm (a partitioning algorithm similar to the k-means algorithm). Microstates are iteratively lumped into macrostates using the k-medoids algorithm, and split again into microstates to iteratively refine the definition of the macrostates.
Voronoi tesselation has also been proposed since it allows a simple construction of macro-states based only on the definition of the cell centers (Vanden-Eijnden and Venturoli [2009a]). See Fig. 3. Given a set of points {x i }, a Voronoi cell V i is defined as: Such cells are convenient since there is a straightforward equation to determine in which cell a given point x is. The centers x i can be obtained in different ways including simulations at high temperature, using nudge elastic bands or the string method, techniques to sample transition pathways, or computing minimum (free) energy pathways. We will show later on that placing centers along minimum energy pathways is often a good choice (see Pan and Roux [2008] for a related discussion). We note that in Schultheis et al. [2005] this issue is altogether avoided by constructing a density-oriented discretization (an attempt to circumvent the curse of dimensionality) based on approximating the probability density in conformational space using a mixture of univariate normal distributions.
Once the macro-states have been defined, one calculates the transition matrix: where τ is the so-called lag-time. The eigenvalues and eigenvectors of this matrix allow computing the different rates in the system, along with mean passage times (Swope et al. [2004a], Chodera et al. [2006]). Park and Pande [2006] use the concept of Shannon entropy to measure and identify non-Markovity. Shannon entropy measures the amount of uncertainty associated with a random variable. Non-markovity can be measured by evaluating the change in uncertainty (Shannon entropy) for a Markov variable X n if one prescribes the previous state X n−1 vs. the last two states X n−1 and X n−2 . Markov chains are such that If the previous state is known: H(X n |X n−1 ) = − x,y P(X n = x, X n−1 = y) ln P(X n = x | X n−1 = y) If the last two states are known: H(X n |X n−1 , X n−2 ) = (39) − x,y,z P(X n = x, X n−1 = y, X n−2 = z) ln P(X n = x | X n−1 = y, X n−2 = z) Measure of non-Markovity: R = H(X n |X n−1 ) − H(X n |X n−1 , X n−2 ) H(X n |X n−1 ) From this definition, R = 0 for Markov chains 2 and R = 1 if the knowledge of X n−2 and X n−1 completely determines X n . This procedure can also be used to refine the definition of macro-states. Another approach to identify and correct non-Markovity is proposed in Nerukh et al. [2010] using the concepts of "computational mechanics" (a term coined by Crutchfield [1994]) and -machines (Crutchfield and Young [1989], Shalizi and Crutchfield [2001]). Huang et al. [2009] discusses the use of generalized ensemble algorithms, e.g., the replica exchange method, parallel tempering or simulated tempering (Hansmann and Okamoto [1999], Sugita [1999], Lyubartsev et al. [1992], Marinari [1992]) and how they can be combined with Markov State models to provide an efficient conformational sampling algorithm. See Bowman [2009] for a discussion of similar ideas and how they have been implemented in the piece of software MSMBUILDER. Chiang et al. [2010] proposes to construct Markov models with hidden states as a way to construct more accurate models with fewer states (hidden Markov models). In such models, the hidden Markov states do not represent individual protein conformations but rather overlapping probabilistic distributions over the conformational space.
Applications of this approach are discussed in several papers including: • a polyphenylacetylene (pPA) 12-mer in explicit solvent for four common organic and aqueous solvents (acetonitrile, chloroform, methanol, and water): Elmer et al. [2005a,b].
• the villin headpiece: Jayachandran et al. [2006].  discusses the application of the software MSMBUILDER to the villin headpiece.
• 4 Non-equilibrium umbrella sampling and reactive trajectory sampling 4.1 Non-equilibrium umbrella sampling In recent years, some approaches have in some sense tried to combine ideas from the previous sections, for example by calculating reactive trajectories or transition pathways from A to B (thereby being exact, contrary to Markov state models whose accuracy depends on the Markov assumption), while using a decomposition of the conformational space Ω into cells. These methods therefore combine the advantages of transition path sampling with the sampling efficiency of MSM, which requires only "local" sampling. The first approach we will discuss is the one proposed by Dinner and co-workers, Warmflash et al. [2007], Dickson et al. [2009b]. The method originates from the method of umbrella sampling, Torrie and Valleau [1977], in the sense that it tries to enhance sampling in poorly sampled region. The latter was broadly speaking adapted to allow modeling non-equilibrium systems. In this case, when computing a reaction rate from A to B, particles are removed from the system each time they reach B and are reinjected in A, thereby creating a steady-state but out of equilibrium situation.
This technique constructs two staggered lattices (using boxes, Warmflash et al. [2007], or following minimum energy pathways, Dickson et al. [2009b]). Simulations are run inside each macro-state (a box in Warmflash et al. [2007]). Then one records when the system attempts to leave a box. At that point, two strategies are applied.
First, we keep track of the number of particles going from box i to j. Each time a particle attempts to go from i to j, the weight of box i, W i , and j, W j are adjusted according to: where T i is the time elapsed in region i (to account for situations in which longer simulations are run in some boxes), T * is some arbitrary time scale to make the equation dimensionally correct, and s is a small parameter used to adjust the rate at which the weights W i s vary. With this equation, the weight of each box converges to its correct steady-state value. Second, one needs to determine which point should be used to reinsert a walker that left box i back into box i. For this walkers that leave another box j and attempt to enter i are saved and, among those, one is picked according to the probability rule where p a is the probability of choosing state a in box i, N ba is the number of crossings from state b in box j(b) = i to state a.
In Warmflash et al. [2007], page 154112-4, end of section G, it is argued that two lattices are needed. Simulations are run in both but walkers that re-enter box i in lattice 1 are chosen with Eq. (43) using data from lattice 2. Otherwise it is claimed that the method is unstable and convergence may not be achieved. The argument put forward is as follows: Suppose, for example, that the weight of a box (B) fluctuates upward. By Eqs. (3) and (4) [Eq. (43) in this manuscript], walkers in neighboring boxes will then be reset to boundary states accessible from B more often. However, if transitions from those states to ones in B are allowed, with some probability, the reset walkers will immediately attempt to enter B and increase its weight further according to Eq. (5) [Eq. (42) in this manuscript]. This positive feedback loop causes the single-lattice scheme to be unstable in simulations to obtain the steady-state probability distribution as a function of multiple variables. The use of two lattices enables boundary states on one lattice to be chosen using the fluxes from the other lattice, which breaks the feedback loop and enables convergence.
This argument is not so clear unfortunately. As a particle leaves i to enter j, the weight W i is reduced by −sW i T * /T i . If the particle reenters i, the weight is increased again by sW j T * /T j . In general if W i has an upward fluctuation, the net result is a reduction in W i . In some cases, walkers do not re-enter i and move to some other box, further reducing W i . As a result, in the scenario mentioned above of a temporary fluctuation upward of the weight, on average, Eq. (42) will slowly reduce the weight of box i and return it to its correct steady-state value.
Following Warmflash et al. [2007] and Dickson et al. [2009b], Vanden-Eijnden and Venturoli [2009b] developed a similar approach but that uses a single partition of space, based on Voronoi cells. The approach is similar with the following differences: • Weights are adjusted based on fluxes between cells. These fluxes, N ij /T i (where N ij is the number of crossings from i to j), are used to solve a linear system that provides an approximation to the steady-state weights W i : Another "global" scheme to adjust the weights (by contrast with the local scheme (42)) is given in Dickson et al. [2009a].
• To pick a re-entry point, one of the boundaries α of cell i is randomly picked using a probability law obtained from the cell flux and steady-state probabilities: where boundary α is the boundary between cell i and j. Although the implementation is different, this is similar to Warmflash et al. [2007].
Based on this, it appears that using two lattices is not necessary and that the scheme correctly works with a single lattice. Also the idea of lattice is no longer discussed in a more recent publication, Dickson and Dinner [2010]. See also Dickson et al. [2011] for an application of this method to unfolding and refolding of RNA. In this paper as well, a single lattice is used. Dickson and Dinner [2010] present some theoretical results regarding non-equilibrium umbrella sampling (an analysis of the convergence of the weights using the local scheme), a comparison with and discussion of forward flux sampling, and recent applications of these methods.

Reactive trajectory sampling
The second method, which is related in some fashion to the previous class of techniques, can be attributed to Huber and Kim [1996]. In this reference, the method is developed assuming that an approximate reaction coordinate has been chosen. However, it is not difficult to extend this approach to a general decomposition of the conformational space Ω in a manner similar to, for example, Vanden-Eijnden and Venturoli [2009a] with Voronoi cells. This method will be discussed in more details below. It consists in running a large number of simulations (or "walkers") in parallel in such a way that a given number of walkers are maintained in each cell or macro-state. Macro-states that are near an energy barrier will tend to be depleted and therefore a strategy is applied to duplicate walkers in this macro-state, in a statistically correct way. This is done by assigning statistical weights to each walker. For example a walker with weight w can be split into two walkers, starting from the same location in Ω, with weights w/2. Conversely, macro-states that are at low energy will tend to become overcrowded and walkers are then removed. If for example we have two walkers with weights w 1 and w 2 , we randomly select one with probabilities (w 1 /(w 1 + w 2 ), w 2 /(w 1 + w 2 )) and assign to it the weight w 1 + w 2 . This approach ensures an efficient sampling of phase space.
In order to calculate a reaction rate, the macro-state corresponding to region B is transformed into a cemetery state, that is any walker that enters this macro-state is removed from the simulation and reinserted in region A. In this fashion, although the simulation is effectively out of equilibrium, the population of walkers is kept constant. This method allows computing all the relevant quantities of interest, such as reaction rates, free energy, metastable states, etc. We note that contrary to Markov state models, this approach does not suffer from non-Markovity errors and that in the limit of infinite sampling it provides an exact answer.
In Zhang et al. [2007], this technique was applied to explore the transition paths ensemble in a unitedresidue model of calmodulin. See also Zhang et al. [2009], . In , it is shown that the method initially developed in Huber and Kim [1996] is really applicable to a much wider class of problems and proposes some generalizations of this procedure.
We mention that a similar technique has been applied to simulated annealing to find minima of rough (or even fractal) functions (see Huber and McCammon [1997]).
Detailed discussion of reactive trajectory sampling. We now discuss in more details the method of Huber and Kim [1996], Zhang et al. [2007 which we rename reactive trajectory sampling method (RTS), in the broader context of macro-state models (e.g., Voronoi cell partitioning). In this approach, systematic errors arising from non-Markovian effects are avoided by directly calculating reactive trajectories from A to B and obtaining the probability flux entering B (or A for the backward rate), Metzner et al. [2006]. When the energy barrier is high, this can be very inefficient since very few trajectories (if any) will make it to B when started from A. However a simple trick allows improving the efficiency of the calculation to the extent that the decay of the statistical errors becomes essentially independent of the energy barrier height.
As before we split the space of possible configurations into cells. Then a large number of random "walkers" are initialized and advanced forward in time. The basic idea is to use a strategy whereby, in cells that get overcrowded (too many walkers), we merge walkers, thereby reducing their numbers, while in cells that are depleted (near transition regions), we split walkers to increase their number. The end goal is to maintain a given target number of walkers in each cell. With such an approach we are able to observe a constant stream of walkers going from A to B (and vice versa) irrespective of the height of the energy barrier. We now explain the details of the method.
Assume we have n w walkers whose position gets updated at each time step. It is possible to resample from these walkers without introducing any bias in the calculation using the following procedure. Each walker, whose position is denoted x i , is assigned a probabilistic weight w i , for example initially equal to 1/n w . A walker can be split into p walkers with weight w i /p. After the split, each walker can be advanced independently. Averages can then be computed using: This equation is always true irrespective of how many times the splitting procedure is applied, or how many steps are performed, as long as the initial position of the walkers is drawn from the equilibrium distribution. This is proved from the fact that the equilibrium distribution is by definition invariant under the dynamics under consideration for x(t).
The reverse operation is possible. Assume we have a group of walkers with weights w 1 , . . . , w nw . Suppose we randomly pick a walker with probabilities w 1 / k w k , . . . , w nw / k w k , and assign to it a weight of i w i . Since the average weight of walker j is ( i w i ) (w j / k w k ), Eq. (46) remains true. This procedure can be used to reduce the number of walkers in a cell.
The algorithm below, called resample, explains the procedure for resampling walkers. Before the procedure resample is called we have several walkers in each cell with varying weights. The procedure resample loops over the cells and select walkers in a way such that the new set of walkers all have the same weight, equal to the average weight of walkers in the cell. This is an important aspect of the method as assigning a constant weight can be proved to be optimal in terms of minimizing the variance and therefore statistical errors. This approach is different from Huber and Kim [1996], which leads to walkers with varying weights and results in a somewhat larger variance. The proof is given below.
1 # This a l g o r i t h m , c a l l e d resample , i s w r i t t e n i n t h e programming l a n g u a g e # python . I t was w r i t t e n by E r i c Darve and J e s u s A. I z a g u i r r e ( U n i v e r s i t y 3 # o f Notre−Dame ) . # I n p u t : l i s t o f w a l k e r s ( l i s t 0 ) and l i s t o f w a l k e r w e i g h t s ( w e i g h t s ) . 5 # I n p u t : t a r g e t number o f w a l k e r s ( n t a r g e t w a l k e r s ) . # Output : l i s t o f w a l k e r s ( l i s t 1 ) and t h e i r w e i g h t s ( n e w w e i g h t s ) .

e i g h t s [ x ] = Wxy
We note that this algorithm terminates since when the last walker is removed from list0 we exactly have W = tw (target weight) so that the while loop does terminate. The maximum number of iterations in the while loop is bounded by n w + tn: the number of walkers in the cell before the procedure starts + the target number of walkers. (Lines 26 through 40 are executed at most tn times, while the lines 42 to 48 are executed at most n w times.) In addition it is apparent from line 33 that the walkers have the same weight at the end. The sorting of the weights on line 13 helps reduce data correlation. Indeed when a walker is split, samples become correlated for some time. The initial sorting makes sure that only walkers whose weight is greater than tw are split. As soon as we have processed all the walkers with weight greater than tw, r stays equal to 1. The reinsertion on line 35 is required to ensure a constant total weight W. The weight that is used, W−r*tw, ensures that the resampling is correct and that the total weight in each cell is unchanged by resample.
Proof of optimality. To simplify the discussion we assume that we have n w particles with weights w i such that i w i = 1. We evolve the system in time such that the walker locations {x i } become uncorrelated. As a consequence, the weights w i are independent from the positions {x i } nw i=1 . We also assume that the walkers do not have any particular order so that the statistics of x i and w i are the same as x j and w j , j = i. The estimator of a particular quantity f isf which tells us thatf is indeed an unbiased estimator of f . The statistical error can be estimated from the variance off : since w 2 = 1/n 2 w . In our algorithm the weights are kept constant (if tn in Algorithm resample is constant) so that Var(w) = 0. The statistical error is therefore minimized. In Huber and Kim [1996], the weights are not constant resulting in a larger statistical error.
An important difference with the Markov model based on cells (coarse states) is that MSM must be run by construction with a known lag time τ . Post-processing is then required to determine whether τ is large enough. After this, the simulation may need to be run again with a larger τ if it is found that memory effects are important. On the contrary, the convergence of RTS is easier to monitor. One simply needs to record the particles that reach B and calculate the average flux. If the error is found to be too large, one can simply continue the simulation to accumulate more statistics, without losing the data already gathered.
In Huber and Kim [1996], walkers that enter B are re-inserted in A thereby ensuring a steady-state system. It is possible to use a slightly different procedure where colors are given to walkers such that when a walker last entered A its color is blue, while it is red if it last entered B. In effect the color changes from blue to red the first time the walker enters B (similarly with A). The population for both colors is kept constant in each macro-state according to the algorithm resample. This algorithm allows computing both the forward and backward rates, the free energy, and the equilibrium distribution of particles (by considering all particles, of any color).

Optimal cells
Even though the method was shown to always converge to the correct answer, the rate of convergence, which depends on the rate at which particles transition from A to B, depends on the choice of cells. We provide some guidelines to help make a good choice of cells. We note that, as explained previously, RTS is always unbiased and is exact in the absence of statistical errors. This is in contrast with milestoning which is exact only when the milestones are iso-surfaces of the committor function. In that case the milestones are called optimal since they minimize systematic errors. Here, RTS is unbiased. The optimal milestones in this context are the ones that minimize the statistical error. Even though the terminology is the same, the meaning is therefore quite different.
RTS ultimately amounts to sampling reactive trajectories from A to B. Reactive trajectories are defined as trajectories that leave A and reach B without reentering A at any point. These trajectories cluster around the minimum (free) energy paths, Maragliano et al. [2006]. These paths correspond to reactive trajectories going from A to B with maximum likelihood, that is the probability density associated with this trajectory is maximum.
To discuss properties of minimum free energy pathways (MFEP), it is convenient to use generalized coordinates (ξ 1 , . . . , ξ p ) to describe the system and use the free energy A(ξ 1 , . . . , ξ p ). This is practically a more useful description and it removes degeneracies such as translation and rotation invariance. For example these generalized coordinates can be chosen as a set of internal coordinates describing the shape or structure of a molecule. If some information is available about the system we can reduce the number of such variables to focus on the variables of interest for the reaction at hand.
The generalized coordinates can have any units, for exampleÅ or deg. This indicates that some kind of non-dimensionalization procedure is required to work with ξ. This non-dimensionalization can be derived in different ways. We shortly describe how this can be done. If one assumes for example a Brownian model for these variables: where R is such that R R T = 2 D, β = (kT ) −1 , and η(t) is a random term with a normal distribution and variance 1, and A is the free energy. The tensor D is the diffusion tensor. Under simplifying assumptions this tensor can be approximated by: where m k is the mass of atom k; τ D is a time scale associated with the rate of decay of the auto-correlation function for dξ/dt. The tensor M ξ is non constant. However we will assume that it can be approximated by its statistical average and that its fluctuations can be neglected. From Eq. (53), we can conclude that along the MFEP we must have that the tangent dξ/ds (where s is some parameterization such as the arc length) is parallel to D∇A or equivalently: This result is somewhat counter-intuitive since we would expect dξ/ds ∝ ∇A but is a result of the metric associated with ξ and defined by M ξ . This suggests normalizing ξ using: With the variablesξ we have the expected relation, along the MFEP: where∇ involves derivatives with respect toξ. Similarly when defining Voronoi cells, the correct distance to use should be consistent with the rate of diffusion and therefore the following distance must be used: We will now assume that we are usingξ instead of ξ, but we will keep the notation ξ for simplicity. Returning to the issue of reactive trajectories and optimal cells, we note that in each cell, walkers tend to accumulate in low energy regions. Consequently we can expect the method to be efficient (practically the statistical errors are small) whenever the regions around the MFEPs are well sampled, that is the low energy region in a cell should overlap as much as possible with the MFEP.
Consider a cell and assume that A is not singular in that cell. This implies that A is minimum on one of the boundaries, which we denote S. Let us assume that the MFEP crosses S at ξ 0 . In order to minimize the statistical errors, we impose the condition that ξ 0 is the point with the highest probability density in the cell. In that case, the boundary must be orthogonal to ∇A at that point, that is the tangent to the MFEP should be orthogonal to the cell boundary. As a note, we point out that along the MFEP, the gradient of On the left panel, the boundary of the cell is normal to the MFEP, ensuring a good sampling along the MFEP. This is illustrated by the most populated region containing ξ 0 . On the right panel, the orientation of the boundary was changed so that the most populated region is now shifted and no longer overlaps with the MFEP. In that case statistical errors are expected to be somewhat larger. In the extreme case where the boundary is parallel with ∇A the simulation converges very slowly.
the committor function is also parallel to ∇A, which implies that locally the boundary S is an iso-surface of the committor function. See Fig. 4. As pointed out previously the method always converges to the correct answer. However the statistical errors are expected to increase as the cells degrade, that is the cell boundaries are no longer orthogonal to ∇A. Then the walkers in a cell start accumulating in regions that are far away from the MFEP. As a result the probability of seeing a walker reach B becomes correspondingly smaller, leading to larger statistical errors. Instead of observing a steady (with small fluctuations) stream of particles with small weights reaching B, we see a more sporadic stream of particles with larger weights. See Fig. 4.

Metastability, convergence, and the multi-colored algorithm
In RTS, walkers initially start from region A and are moved out of the simulation when they reach region B. If nothing is done, the total weight slowly diminishes. When the reaction rate is slow, the simulation remains accurate. However, when multiple rates are comparable, this may lead to biased results. One option to prevent the total weight from diminishing is to reinsert a walker in basin A after it enters basin B. Another possible strategy, as discussed earlier, with similar efficiency but somewhat easier to implement, is to use walkers with two different colors, for example blue and red. By definition, red (resp. blue) walkers are those that have last visited A (resp. B). Whenever a blue walker enters B, its color changes to blue, and vice versa. The resampling algorithm is applied to each color separately so that each bin contains the same amount of walkers of each color. This approach allows computing the forward and backward rates, and the free energy.
The efficiency of RTS degrades when there are other slow reactions rates (values of λ k that are small) in the system, that is other metastable regions inĀ orB, Noé et al. [2007]. In that case, the convergence of the algorithm is limited by the rate at which walkers equilibrate in regionĀ andB, which is given by λ −1 3 . This leads to long correlation times for the measured fluxes, and therefore overall slow decay of the statistical errors. This can be addressed by assigning appropriate weights for the walkers at t = 0. One possibility is to start running RTS for a while and calculate fluxes between cells, given by the weights of walkers moving from cell i to j after a time step divided by the total weight of walkers in cell i. From the flux matrix, the steady-state weights of cells can be computed. These weights can then be used to adjust the weights of walkers in each cell. In principle this needs to be iterated until convergence. However only one or two iterations are typically needed. The remaining deviations from steady-state can be later on reduced by the production RTS run. The goal of this procedure is simply to improve the initial weights in each cell in order to bypass the initial slow convergence.
The key property of this procedure is that the accuracy with which fluxes can be computed is independent of the presence of metastability and depends primarily on the sampling inside each cell. Overall convergence is a global property, which is why a few iterations may be required, but since this is followed by the production RTS run, this part of the calculation only requires low accuracy and few iterations. This approach is similar to a method described in , called "Enhanced weighted ensemble attainment of steady state." Another approach can be applied to address this shortcoming. It consists in using more than 2 colors. To explain this in a simple fashion we return to the milestoning framework of Section 2.6. RTS can be reinterpreted in terms of the milestoning framework. The difficulty in the optimal milestoning approach described in Section 2.6 is that the cells must be such that their boundaries are iso-surfaces of the committor function. Even though this is in general difficult to realize, there is a case for which the problem is simplified. Consider the following three cells: cell V 0 enclosing A, cell V 1 enclosing B. Let us assume that we have another metastable basin C. The same reasoning can be extended to an arbitrary number of basins. We define V 2 as a cell enclosing C and V 3 as the complement: We define S 0 as the interface between V 0 and V 3 , and similarly for S 1 , and S 2 . In the multi-coloring algorithm, each cell is assigned a color, say 0 is blue, 1 is red, and 2 is green. Each time a blue particle enters cell V 2 , its color changes to green, and similarly for the other colors. The total weight of all walkers is therefore constant. We use the resampling algorithm to maintain a constant number of walkers in each cell, separately for each color. We then measure the mean flux of particles turning blue to green, etc.
Eq. (27) still holds for this approach. The fluxes F ij , from cell i to j, are simply related to the probability matrix P through: F = ∆t −1 P . Eq. (29) is not needed for this approach. We can directly obtain P from the flux values F , computed using RTS.
In general the milestone S 2 is not an iso-surface of the committor. However, following the proof that optimal milestoning gives exact rates, we will have proved that the rate with multi-coloring is exact if we show that: p ij , the probability to cross j after i, and τ i , the mean time before crossing another milestone, are independent of the previous milestone k that was crossed (see page 14 for the proof in the optimal milestoning case). This property is in fact true for the multi-coloring approach because of our choice for S 2 and the fact that the equilibration time for V 2 is very small. The key assumption is that V 2 must be associated with a metastable state (minimum energy basin) so that the relaxation time in V 2 is small compared to the mean escape time.
This shows that the mean passage time predicted using Eq. (27) in the multi-coloring framework is for all practical purposes exact, with no systematic bias. This approach allows considering the case of multiple slow rates, with no significant degradation in efficiency. The computational cost merely grows with the number of colors, but is independent of the degree of metastability (the values of the first λ k , k = 1, 2, etc). This approach has some conceptual similarities with the technique of core sets of . Their analysis of accuracy (section E, pp. 204105-7) carries over to RTS with multiple colors. In addition, if ones applies the Galerkin discretization approach from  to RTS with m + 1 colors, one can calculate the slowest m rates in the system (or phenomenological rates, see  pp. 204105-4, section F) which may be of interest for certain applications.
As a final note, we point out that the method is embarrassingly parallel and can be easily implemented on a parallel machine. This is important as this allows making only few changes to a serial (sequential) molecular dynamics code to make it run efficiently on a parallel cluster, with RTS, without having to parallelize the core of the code.

Analysis of statistical errors in Markov state models
We now discuss some mathematical results for Markov state models. In particular we will analyze the sensitivity of the eigenvalues to perturbations in the transition matrix. This analysis will lead to estimates for the statistical errors in the method. This will also lead to an analysis of the systematic errors, due to the finite lag time τ (the length of the short trajectories used to build the Markov state model).
Several papers have discussed error analysis in the context of conformation dynamics (Sarich et al. [2010], Prinz et al. [2011]) or Markov state models (Hinrichs and Pande [2007], Singhal and Pande [2005], Hinrichs [2007]). Metzner et al. [2009a] takes a different approach to error analysis by considering a method that generates random transition matrices and as a result can estimate errors in various quantities computed from the transition matrix. In Metzner et al. [2009a], it is argued that this approach is more accurate since it does not rely on Taylor expansions to approximate the impact of small variations of the numerical transition matrix from the exact matrix. In that sense this is a more direct estimate of the statistical errors and how they impact various quantities of interest (stationary distribution, eigenvalues, committor function, etc). Darve et al. [2009] discusses the effect of memory in building coarse grained models in the context of the Mori-Zwanzig formalism. Although many results presented in this paper can be extended to more general stochastic equations, we are going to focus on the relatively simpler case of Brownian dynamics (Ermak and McCammon [1978], p. 1355): where W (t) is a Wiener process (see p. 66 in Gardiner [1997]), D B is the diffusion tensor, β −1 = k B T , and

Eigenvectors and eigenvalues of the transition matrix
The rate is typically computed by considering the eigenvalues of the transition matrix: where we assume that the dynamics is given by (60). This matrix is basically used to construct a Markov state model of the system. From the stochastic equation (60), we can define the conditional probability ρ(x, t|x 0 , 0), which is the probability to be at x at time t if the system was at x 0 at time 0. This probability can be expanded in terms of the eigenfunctions ρ k (x) of the forward Fokker-Planck equation (for Eq. (60)), and the eigenfunctions ψ k (x 0 ) of the backward Fokker-Planck equation (Gardiner [1997], p. 165): where λ k are real and positive eigenvalues (the two sets of eigenfunctions are associated with the same eigenvalues). See p. 32 in Schütte [1999], p. 174105-5 in Prinz et al. [2011], and p. 166 in Gardiner [1997]. Since e −λ k t is an eigenvalue of the kernel ρ(x, t|x 0 , 0): We will denote: where ρ(x) is the equilibrium distribution of the system. We sort the λ k in increasing order. For most systems, there is a single eigenvalue exp(−τ λ 1 ) equal to 1 (λ 1 = 0) and the corresponding eigenvector is the stationary distribution ρ(x). We are interested in estimating λ 2 by computing the eigenvalues µ 2 of the matrix P ij , and using λ k ∼ − ln(µ k )/τ . In general, λ 2 and − ln(µ 2 )/τ differ leading to inaccurate estimates. However under certain assumptions, which will be discussed, − ln(µ 2 )/τ provides an accurate estimate. In these circumstances, the Markov assumption made in building the model becomes accurate. The term P ij can then be written as: In the rest of this paper we sometimes have to make a distinction between the two minimum energy regions A and B and a partitioning of the space Ω into two metastable regionsĀ andB (see Fig. 5). Typically A is defined as a small region around a stable conformation of interest (the reactant state), and similar for B (the product state). The regionsĀ andB are defined as metastable regions, that is the rate of transition between these regions is the smallest among all other pairs of sets.
Properties of the eigenvectors. The second eigenvector ρ 2 (x) can be used to define a partition into two regionsĀ andB, with A ⊂Ā, B ⊂B. We will provide a more rigorous analysis later on but roughly speaking, the function ρ 2 /ρ is nearly constant in two regions, which defineĀ andB. This function changes sign between these two regions. The narrow region where ρ 2 /ρ ∼ 0 defines the transition region between these two metastable regions. See for example Chapter 6, p. 91-119 in Hill and Moore [2004]. The molecule is assumed to have a high probability of being in region A when inĀ, and similarly for B.
We now discuss in more details the properties of the eigenvectors. If we assume that the relaxation time in basinĀ andB is short compared to the reaction rate, i.e., 0 < λ 2 λ k for k > 2, then for x 0 ∈Ā, x ∈B, and λ −1 from which we see that ρ 2 ∝ ρ inB (and similarlyĀ), and ψ 2 must be approximately constant inĀ (and similarlyB). Moreover, we have ρ 2 (x)dx = 0 from which (with the appropriate normalization): Finally we have the general relation A short proof is provided in the appendix (see Proof 1). The eigenvectors ρ 2 and ψ 2 are depicted in Fig. 6. Near the transition region ψ 2 (x) ∼ constant no longer holds, and in fact, the function ψ 2 varies rapidly near the transition region. The committor function π(x) is defined as the probability to reach region B x ψ 2 = ρ 2 /ρ ρ(x) (ρ(B)/ρ(Ā)) 1/2 Midpoint −(ρ(Ā)/ρ(B)) 1/2 Transition point π(x) = 1/2 Points go to setĀ with probability ρ(Ā) RegionĀ RegionB Figure 6: Definition of the metastable regionsĀ andB and the second eigenvector ρ 2 .
before reaching region A starting from x. By definition, at the transition region, π(x) = 1/2. The function ψ 2 has a natural interpretation in terms of a committor function (see Ren et al. [2005]) and See the appendix for a proof of this result (Proof 2). Therefore, from Eq. (71), the transition region corresponds to ψ 2 (x) = 1/2 (ψ 2 (A)+ψ 2 (B)), where ψ 2 (A) denotes the constant value of ψ 2 in A (similarly for B). Therefore, at the transition point: The transition point is shown on Fig. 6. The hypersurface corresponding to ρ 2 (x) = 0 is located nearby but corresponds to a slightly different situation. Points on this hypersurface are not at the transition region but rather are such that they move to basinĀ with probability ρ(Ā) and toB with probability ρ(B). (This is true only after a short time τ , with τ λ −1 k , k > 2.) This is shown on Fig. 6.
Forward and backward rates. Now that we have derived estimates for ρ 2 (x), we can clarify the relation between different rates. To calculate the forward rate, we construct a linear combination ρ f (x, 0) of ρ and ρ 2 such that: With these conditions we find that the unique solution is: Since ρ and ρ 2 are eigenvectors: Therefore the forward rate from A to B is equal to λ 2 ρ(B). Similarly the rate from B to A is equal to λ 2 ρ(Ā): rate A→B = λ 2 ρ(B), rate B→A = λ 2 ρ(Ā).

Sensitivity of eigenvalues and systematic errors
Sensitivity of eigenvalue. We derive a general result regarding the sensitivity of an eigenvalue to perturbations in the matrix entries. We consider a matrix P and assume that we have a full set of eigenvalues. The eigenvalue of interest is µ 2 and: and Λ(P ) is a diagonal matrix containing the eigenvalues. We denoteP the exact matrix P and consider small deviations ∆P = P −P . We assume that in some neighborhood aroundP the matrices R(P ), S(P ), and Λ(P ) are differentiable. Calculate the derivative with respect to one of the entries P ij : We denote e i a row vector such that [e i ] j = δ ij , r 2 = R 2,: (second row of R), s 2 = S 2,: . If we evaluate the partial derivative at P =P , and multiply Eq. (78) to the left by s 2 and to the right by r T 2 we get: because s 2 R T (Λ − µ 2 I) = 0 and (Λ − µ 2 I) S r T 2 = 0. This matrix will be plotted later on, in Fig. 9.
Systematic error due to the lag time τ . We will use this result regarding the sensitivity of µ 2 to study the systematic error or bias using the Markov state model. To simplify the discussion, we will define a fine matrix, which is assumed to accurately capture the continuous dynamics using the Brownian dynamics (60): over some fine states. We consider that P f gives the exact rate (this can be achieved using states that are fine enough). We could work with ρ(x, t|x 0 , 0) directly but it is easier to discuss the results using P f . The matrix P can be written in terms of P f : From P , we can reconstruct an approximation P c of P f where an entry in P is mapped to a block in P c : with k ∈ V i , l ∈ V j , and where |V j | is the number of fine states in V j . We have Π I = I, the identity operator. We will use our result for the sensitivity of µ 2 with ∆P = P c − P f . We now show that every eigenvalue of P is an eigenvalue of P c , which allows us to study the eigenvalues of P c instead of P . Let us defines 2 the left eigenvector of P with eigenvalueμ 2 . Defines c 2 : Therefores c 2 is a left eigenvector of P c with eigenvalueμ 2 . Using Eq. (79), the sensitivity of the second eigenvalue is: where s f 2 is the left eigenvector, and r f 2 is the right eigenvector of P f . We are going to use the following linear approximation:μ We start by studying the systematic error at long lag times. Then: In the following we will make the following approximations: • ρ(x) is negligibly small near the transition region.
• ψ 2 (x) is nearly constant away from the transition region.
In practice this is not true but the error due to these approximations is typically much smaller than O(e −λ3τ ) and so we will simply ignore it. 3 As a consequence, from Eq. (88) In the long lag time assumption, the choice of cells is not important. However, there is one property, which must be satisfied which is that no cell V j can overlap significantly with bothĀ andB, formally: either Ā ∩Vj ρ(x)dx or B ∩Vj ρ(x)dx must be negligible. In that case we have the following three possibilities, assuming that k is away from the transition region: if V j has support in a region where ρ is negligible (transition region).
with an error of order O(e −λ3τ ). This results from the fact that by construction P c kl is constant inside l ∈ V j and that l∈Vj P c kl = l∈Vj P f kl + O(e −λ3τ ). Therefore: At long lag times the error in the eigenvalue is therefore: At short lag times, we have a small systematic error provided the cells are chosen adequately. Let us assume that inside each cell [r f 2 ] l is nearly constant. 4 We will therefore denote [r f 2 ] j the value of [r f 2 ] l for l ∈ V j in the equation below. The fluctuations of [r f 2 ] l in a given cell are assumed to be of order ε. We also have [s (70)], so that: from the definition of P c and P . This suggests the following choice of cell V i that satisfies our assumption: The interpretation is therefore that the cells provide a fine partitioning based on the iso-surfaces of the eigenvector ψ 2 (x) or equivalently the committor function π, and we can equivalently write: In that case:μ This result is consistent with Eq. (35) in Prinz et al. [2011] which gives a bound on the error that depends on the eigenfunction approximation error when projecting onto the cells. This requirement is very important as constructing cells with small volume is prohibitive in high dimension. However the definition (96) is much less restrictive in terms of computational cost since it requires refining only along a single direction, given by ∇ψ 2 . The dimensionality of the problem has little impact on the number of cells that are required for an accurate calculation.
The main caveat in this discussion is that computing or even approximating ψ 2 or the committor function π is very challenging, and consequently this approach remains difficult to use in practice.

Statistical errors
Even though we have outlined methods to reduce the Markovian approximation error, e.g., by refining the spatial discretization, it remains the case in practice that cells need to have a relatively large volume and the committor function is difficult to approximate at best. As a result, a long lag time is more or less necessary for an accurate estimate. We will now discuss how statistical errors vary as a function of the lag time.
At small lag times, the statistics are typically expected to be reasonably accurate and the statistical error can be made satisfactorily small. At long lag times however the situation worsens. In this section, we will again use Eq. (79) that expresses the sensitivity of the eigenvalue µ 2 to perturbations in the matrix entries P ij , this time due to statistical errors. For this analysis we will assume that the macro states are fine enough, or that τ is large enough, so that P (2τ ) = P (τ ) 2 . This simplifies the analysis, but the results mostly carry over to the case where non-Markovian effects are important.
The sensitivity of µ 2 with respect to the matrix entries can be combined with the statistical errors in the entries P ij ({P i,1 , . . . , P i,n cell } is a multinomial distribution, Hinrichs and Pande [2007], Singhal and Pande [2005]) to obtain an expression for the statistical error in the form: where σ i (r 2 ) is the standard deviation of r 2 computed using the probabilities P i,: (row i of P ). The rate is given by λ 2 = − ln(µ 2 )/τ and the relative error can be estimated using σ(λ 2 )/λ 2 : We may be interested in minimizing the statistical error, given a computational cost. The cost is proportional to nτ /∆t, the number of samples multiplied by the length of the trajectories with lag-time τ . We express the error in the form: With our assumption that the states are fine enough, we have that [s 2 ] i ≈ ρ i,2 and [r 2 ] i ≈ ψ 2 i . Therefore if we assume that λ 2 λ 3 , we have that r 2 is nearly constant. Hence σ 2 i (r 2 ) is expected to be small and σ(λ 2 )/λ 2 can remain bounded even as λ 2 → 0.
Again, using the assumption that the states are fine enough, the eigenvectors s 2 and r 2 of P are independent of τ . At short times, the probability P i,: is concentrated around the diagonal (see Fig. 9). In fact for τ = 0 the matrix is equal to the identity and σ i (r 2 ) = 0. As τ increases the spread of the entries in P i,: becomes larger and as a consequence σ i (r 2 ) must increase with τ . The factor 1/(e −λ2τ √ τ ) on the denominator in Eq. (103) results in an initial increase of the error at small τ , then a plateau is reached, and as λ 2 τ 1 the error starts increasing again (τ is at this point large compared to the relaxation rate of the system). See Fig. 9.
An interpretation is that as τ λ −1 3 , the system has time to relax within basinĀ orB. Therefore independent of where the system starts from, the states that were started in basinĀ will be distributed as ∼ ρ(V i )/ρ(Ā) while the states inB are distributed as ∼ ρ(V i )/ρ(B). Therefore in that regime, the method degenerates to a direct calculation of the rate where trajectories are initiated in basinĀ (resp.B) and we observe how many transitions to basinB (resp.Ā) occur. For this type of calculation, statistical errors are large whenever λ 2 τ is very small.
More numerical results will be shown later on, but to illustrate the point above we present a simple example. Fig. 7 shows a 1D system with 0 ≤ x ≤ 1. A random walker is moving between discrete states. The probability to attempt a move to the left is 0.25 (same for the right). A Metropolis criterion is used to accept or reject this move (p. 39, and Frenkel and Smit [1996] p. 27) using: The eigenvalues of the matrix are shown on Fig. 8, along with the decay of exp(−λ k τ ) as a function of τ for k = 2, 3, 4. In Fig. 9, the matrix P is shown along with −log10(P 100 ). This shows how the matrix P τ progressively goes from a tri-diagonal form to a rank-2 matrix given by the first two eigenvectors. Recall the assumptions from page 34 regarding ρ(x) and ψ 2 (x). In Fig. 9, this can be seen from the fact that in the top and bottom portions of the matrix the entries in the columns are nearly constant. In the center of the matrix, entries in the column vary rapidly ( ψ 2 f k is changing sign) but this is also the region where ρ is small. This can be seen by observing that the columns in the center have small values (∼ 10 −4 ). See also Eq. (88). The same behavior can be observed in Fig. 7 where for x < 0.4, ψ 2 is constant; for 0.4 ≤ x ≤ 0.6, ρ is small; for x > 0.6, ψ 2 becomes constant again. Fig. 10 shows the relative statistical error calculated as 2σ(λ 2 )/λ 2 . As τ increases the statistical error increases because the system has more time to relax in each basin, thereby reducing the computational benefit of using coarse states.
The overall behavior is therefore as follows. At small lag times, we have non-Markovian effects leading to systematic errors. As τ increases, the terms e −λ k τ , k > 2, become small when λ k τ 1, such that systematic non-Markovian effects progressively disappear. As this happens, the statistical errors increase because the system has more time to relax in each basin (see Fig. 10). When we reach τ λ −1 3 (around τ ∼ 40), we see a plateau. Systematic errors are now negligible but the statistical errors are large.
The implication is that it may be difficult in practice to apply these methods accurately. The reasoning above shows that if the cells systematic errors at small lag times and large statistical errors at large lag times, and a trade-off must be found between these two extremes to maximize the efficiency of the calculation. This issue can be mitigated by choosing cells V i such that ψ 2 (x) ∼ constant in each cell, although this can lead to a large computational cost if the number of cells becomes too large or may be intractable if the committor function cannot be well approximated.
The advantage of RTS is that none of these issues are present. Systematic errors are absent and convergence is easy to monitor. There is no lag time τ that needs to be adjusted to control the accuracy and computational cost. One benefit of MSM though is that it requires only independent sampling inside each macro-state, making the calculation embarrassingly parallel. In contrast, RTS requires a "global" convergence of the macro-state weights. This is a necessary trade-off however resulting from the fact that the Markov approximation is not used and therefore transition pathways, in effect, must be computed.
6 Some numerical benchmarks for Markov state models and reactive trajectory sampling

Set up of numerical simulations
To illustrate the performance and numerical behavior of Markov state models (MSM) and the reactive trajectory sampling method (RTS), we present results from some numerical simulations on simple problems in 1D and 2D. They are meant to illustrate the concepts introduced in the previous sections. These results can not really be used as benchmarks of the various methods presented, as conclusions using actual models from bio-chemistry will inevitably be different. Nonetheless they illustrate the key concepts and confirm to some extent the validity of the analysis. We will consider two types of dynamics. In 1D, we will consider Brownian Dynamics as introduced in Eq. (60). Using a time step ∆t, the numerical integrator is given by Ermak and McCammon [1978]: where D B is the diffusion matrix, W is a vector of independent standard Gaussian random variables (mean µ = 0, standard deviation σ = 1). A special difficulty of this discretization is that depending on the choice of time step, the walker x n+1 may (infrequently) jump over many macro states. In fact there is a finite chance to jump in one step directly from A to B. This is an artifact of the discretization. This case is a little artificial but does highlight a limitation of the method, which is that the statistical error typically increases when such "long-range" jumps are allowed. In practice however, in the context of molecular systems, this is not an issue as the system typically diffuses slowly and can only cross to neighboring cells during a single time step. In 2D, the walkers are moving on a 2D Cartesian grid and we consider the Metropolis algorithm (see Frenkel and Smit [1996] p. 27). The scheme is outlined in algorithm 1.
Algorithm 1: Metropolis Algorithm /* Walkers are moving on a Cartesian grid in dimension d.
*/ 1 while more samples are required do 2 Move walker to an adjacent position x new with uniform probability of 1/2 d (d-dimensional random walk.) 3 Accept this move with probability If move is rejected, stay at x old , otherwise move to x new .
The exact rate is obtained by computing the eigenvalues of Q ij = P(walker at j the next step | walker currently at i) ≈ ρ(x j , ∆t|x i , 0) In 2D the system already has discrete states so Q is well-defined. In 1D, we discretize the interval of interest in order to define the fine states i. The coarse transition probability matrix P represents the transition probability between the coarse cells or macro states. The definition is simply as follows: This matrix can be computed analytically in our examples because everything is low dimensional. For example we can use: where P(y|x) is obtained from Q and P(x|V i ) is known exactly. Now given a measurement of P , one can compute the second eigenvalue µ 2 and let − ln(µ 2 )/τ be an approximation to λ 2 . The RTS algorithm proceeds as follows.
In our examples we only have two basins and therefore only use two colors, red and blue. Each time a red particle enters B, its color changes to blue, and vice versa. We use Algorithm resample to maintain in each cell or macro state a constant number of walkers of each color. At each step n, we calculate J b→r n = ∆t −1 Sum of the weights of all the blue particles who turn red during step n Sum of the weights of all the blue particles at step n The rate from A to B is then given by: The reverse rate, rate B→A , is obtained similarly.

1D Brownian Dynamics Setup
The 1D simulation takes place in the domain Ω = [−10, 10]. The system parameters are β −1 = 0.2, D = 0.06 (diffusivity), ∆x = 0.03 (fine state discretization), ∆t = 0.03, n walkers per cell = 10, and n cell = 32 (coarse states). We define:Ā = [−10, 0], A = [−7, −5],B = [0, 10], and B = [5,7]. The potential U (x) is The exact rates computed using the fine states (with spacing ∆x and time step ∆t) are rate A→B = 1.59 × 10 −8 and rate B→A = 6.70 × 10 −11 . The overall setup is illustrated in figure 11. Figure 11: A snapshot of the 1D simulation setup. The green line shows the potential U (x) and the red and blue lines show the positions and weights of the walkers. The weight has large variations and as a result we used a non-linear scaling to determine the length of the short vertical lines with l = β(log(αw + 1) + 1), where l is the plotted length, w is the weight of the walker, and α, β are constants. (A simple log-scale is not suitable since it would lead to negative vertical line lengths.) The initialization of walkers can impact the convergence time of the rates. A simple choice is to first scatter red walkers uniformly in region A and blue walkers in region B and then assign weight to the walkers proportional to ρ(x). The resulting distribution of red and blue walkers is not the steady-state distribution as the steady-state distribution has red walkers in region B with small weights. However this is a good approximation of this steady-state distribution. The computed rate from the algorithm will not be correct until the distribution of the walkers has converged to steady-state. Fig. 14a shows results using the coarse grained matrix and estimating the rate using the second eigenvalue µ 2 = e −λ2τ . As predicted by the theoretical analysis of the previous sections we observe two trends. At small lag times τ , the rate is over-predicted because of non-Markovian effects. At longer lag times, the statistical error increases. A good value of the lag time in this example is probably around τ ≈ 300, for which the systematic and statistical errors are both reasonably small.

Discussion
The rates obtained using the reactive trajectory sampling method (RTS) is shown on Fig. 12. This case is actually relatively difficult as the system is allowed to jump across multiple cells in one time step. As a result there might be some infrequent jump of particles with large weights from basin A to B (or vice versa). This results in larger statistical errors, which leads to a large sampling to reach an acceptable accuracy. The initial bias at small times (underestimation of the rate) is caused by the initial distribution of particles, which needs to equilibrate before correct statistics can be obtained. In practice a more judicious choice of initial distribution for the particles can reduce this initial equilibration. : (x + 1) 2 + y 2 ≤ 0.4} (see Fig. 13b). The potential in use is U (x, y) = e −x 2 + y 2 and is illustrated in Fig. 13a.
The domain Ω is divided into n cell cells where the boundaries are equally spaced lines with a slant angle of θ. Fig. 13b is a snapshot of the simulation with a value of θ = 20 • where the boundary of the cells are depicted with the black slanted lines. The time stepping was done using the metropolis algorithm.
If basin A was defined as a vertical slab, the iso-surfaces of the committor function would be vertical lines because of the special form of U in which x and y are basically decoupled. Since A and B are defined as half circles, the choice of θ = 0 • makes the cell boundaries close to the iso-surfaces, and thus should be the optimal choice. We will see that as the angle θ increases the statistical error in the prediction increases.
The exact rates are rate A→B = rate B→A = 5.9 × 10 −7 .

Discussion
Results obtained using the coarse grained matrix are shown in Fig. 14b, 14c, 14d. As before the non-Markovian effects are visible. From the theoretical analysis, we established that non-Markovian effects are reduced when the cells boundary approximate the iso-surface of the committor function and the cells are narrow in the direction ∇π. Here the width of the cells is kept fix but we vary the angle of the cells boundaries, which are straight lines. At θ = 0, the lines are relatively good approximation of the iso-surfaces and therefore memory effects are small. As θ increases, the deviation from the iso-surfaces increases and we observe longer lag times. Note that the cell width is the same so that the degradation does not correspond to a coarsening of the cell but rather a poor choice of their geometry. In complex examples from bio-chemistry we expect that the choice of cells is far from ideal since it is difficult to guess the iso-committor surfaces and therefore the situation in Fig. 14c and 14d is somewhat representative. We also note that in addition to larger memory effects, the statistical errors increase with θ.
The rates computed using RTS are shown in Fig. 15a-15e. We see that as expected the rate always converges towards the correct value. As θ increases the statistical errors become larger and at θ = 40 the method is no longer very efficient. Even though it is shown, the choice θ = 90 would basically ruin the  method since the cells would no longer provide any sampling enhancement. In that case, all the cells would have horizontal boundaries, parallel to x. Fig. 13b illustrates the point made in Section 4.3. The minimum energy path in this case is the centerline y = 0. As the cells become tilted particles tend to move away from the centerline. In this case, on the left side, particles tend to move up and on the right side particles tend to move down. Roughly speaking, particles accumulate near the point where ∇A is orthogonal to the cell boundary. As θ increases, this point moves away from y = 0. As a result fewer particles are found near the center lines, even though these are the particles that make the largest contributions to the rate. This translates into infrequent events in which a particle with large weight changes color. At θ = 20 this effect is small but it increases with θ. At θ = 90, we would see very few particles changing color but with a very large weight. The rate would therefore be still correctly calculated, however with a large standard deviation.

Conclusion
This paper has discussed a number of approaches to calculate reaction rates. We reviewed the reactive flux approach, transition path sampling, transition interface sampling, forward flux sampling, conformation dynamics, Markov state models, non-equilibrium umbrella sampling, and an extension of weighted ensemble Brownian dynamics (renamed reactive trajectory sampling in this paper). These methods differ in their assumptions and computational cost.
The reactive flux approach is probably the least computationally expensive but requires a fair amount of knowledge about the system (transition pathways, transition region, location of saddle point, etc). Transition interface sampling (TIS) and forward flux sampling (FFS) both rely on a reaction coordinate ξ or an order parameter (variable that increases monotonically with ξ). In cases where the mechanism is unknown or multiple pathways are present such methods may be more difficult to apply.
Markov state models (MSM) attempt to alleviate some of these problems by considering a general partitioning of space in macro-states. They can therefore more easily accommodate multiple pathways and complex reaction mechanisms. MSM depends on a lag time τ , equal to the length of trajectories that need to be run. At short times τ , significant non-Markovian effects are present, resulting in a systematic bias,  As θ increases the method shows no systematic bias but statistical errors increase. The initial relaxation to the exact rate is a result of the initial distribution of particles. The legend is the same for all plots and is only shown for the first plot.
while at long times, statistical errors increase. MSM are "relatively" easy to set up and calculate as they involve running a large number of independent short trajectories. It is therefore embarrassingly parallel and can be run efficiently even on loosely connected computers, for example with cloud computing. Some methods attempt to improve on transition path sampling methods (TIS, FFS) and MSM. They include the non-equilibrium umbrella sampling and weighted ensemble Brownian dynamics, which we called reactive trajectory sampling (RTS) when extended to the case of general macro-states. In this case the assumption of Markovian dynamics can be relaxed. This comes at the cost of a global convergence of statistics across macro-states. For example in RTS and non-equilibrium umbrella sampling, the weights of the macro-states need to be converged, typically using a fixed-point iteration scheme. Although this is typically fast, this does lead in general to a computational cost that is larger than with MSM. In addition, some amount of global communication is required at regular intervals (e.g., for the resample algorithm in RTS), although the information that needs to be communicated is minimal. It includes the weights of walkers w i and macro-state index. Depending on the setup, particle coordinates (x) may need to be communicated in some cases.
Techniques exist to construct optimal macro-states or improve an initial guess. They rely for example on computing minimum energy pathways (MFEP). Although there have been many successful implementations, this remains a challenging problem. For example, computing all the MFEPs that make significant contributions to the reaction rate and building macro-states from this data remains a challenge.
Nevertheless, these methods offer promising avenues to calculating reaction rates, transition states and reaction mechanisms. They make excellent use of modern parallel computers as most of these methods involve running a large number of fairly independent trajectories (perhaps with a small amount of global communication required). This puts less pressure on developing software that can run long trajectories using many processors. The efficiency of these methods is independent of the degree of metastability of the system (Perron cluster) and their convergence is often dictated by the shorter mixing time scales inside each metastable basin (at least in an optimal set up of the method).
Acknowledgements. The authors gratefully acknowledge the work of Jesús A. Izaguirre and Haoyun "Michelle" Feng (University of Notre-Dame) who made several suggestions to improve the manuscript.

A Technical proofs
These proofs can be found in classical textbooks such as Gardiner [1997]. We provide them here as they can be helpful to understand some of the results and are also sufficiently simple to be succinctly explained.
Proof 1. We prove that ρ k (x) = ψ k (x)ρ(x) Consider simply Eq. (60) in one dimension (although the proof can be extended to the multi-dimensional case): The function ρ(x, t|x 0 , 0) satisfies the forward Chapman-Kolmogorov differential equation: From the eigenfunction expansion we therefore have: Similarly the backward equation is satisfied. In this case we write: ρ(x, t|x 0 , 0) = ρ(x, 0|x 0 , −t) and ∂ρ(x, t|x 0 , 0) This leads to: For the equilibrium density ρ (which is equal to ρ 1 ): Consider now ψ k (x)ρ(x). From Eqns. (113), (116), and (117), we can prove that: so that: Proof 2. We prove that We again consider simply a one-dimensional system. Then, the forward Chapman-Kolmogorov equation can be re-written as: ∂ρ where J is interpreted as a probability flux. Assuming absorbing boundary conditions at A and B, the probability to reach B before A is obtained by integrating the flux from 0 to ∞: assuming that B = [b, ∞). From the definition of the flux: This time, we use the backward Chapman-Kolmogorov equation in order to obtain derivatives with respect to x: Since for x = b, J(b, 0|x, 0) = 0, we have: This is the differential equation satisfied by the committor function. This expression can be extended to the multidimensional setting. The boundary conditions are π(a) = 0 and π(b) = 1. The function ψ 2 satisfies a slightly different equation: However with λ 2 very small we can conclude that ψ 2 satisfies nearly the same differential equation as π.
From the boundary conditions, we conclude that: