Efficient estimation of rare-event kinetics

The efficient calculation of rare-event kinetics in complex dynamical systems, such as the rate and pathways of ligand dissociation from a protein, is a generally unsolved problem. Markov state models can systematically integrate ensembles of short simulations and thus effectively parallelize the computational effort, but the rare events of interest still need to be spontaneously sampled in the data. Enhanced sampling approaches, such as parallel tempering or umbrella sampling, can accelerate the computation of equilibrium expectations massively - but sacrifice the ability to compute dynamical expectations. In this work we establish a principle to combine knowledge of the equilibrium distribution with kinetics from fast"downhill"relaxation trajectories using reversible Markov models. This approach is general as it does not invoke any specific dynamical model, and can provide accurate estimates of the rare event kinetics. Large gains in sampling efficiency can be achieved whenever one direction of the process occurs more rapid than its reverse, making the approach especially attractive for downhill processes such as folding and binding in biomolecules.


I. INTRODUCTION
A wide range of biological or physico-chemical systems exhibit rare-event kinetics, consisting of rare-transitions between a couple of long-lived (meta-stable) states. Examples are protein-folding, protein ligand association, and nucleation processes. Meta-stability can be found in any system in which states of minimum energy are separated by barriers higher than the average thermal energy.
A thorough understanding of such systems encompasses the kinetics of the rare-events, e.g. rates and transition pathways. Obtaining reliable estimates for such systems is notoriously difficult: The simulation time needs to exceed the longest waiting time, resulting in a sampling problem.
In recent years, Markov state models (MSMs) [1][2][3][4][5][6][7][8] and their practical applicability through software [9,10] have become a key technology for computing kinetics of complex rare-event systems. A well-constructed MSM separates the kinetically distinct states and captures their transition rates or probabilities. With a suitable choice of state space discretization and lag-time, kinetics can be approximated with high numerical accuracy [8,11]. MSMs can be straightforwardly interpreted and analyzed with Markov chain theory.
A holistic picture of the mechanisms of rare-event processes is possible using Transition Path Theory [12,13], such as demonstrated for protein folding [14,15] or protein-ligand binding [16,17]. MSMs can somewhat alleviate the sampling problem by virtue of the fact that they can be estimated from short simulations produced in parallel [14,15,18], thus avoiding the need for single long trajectories [19]. However, the rare-events of interest must be sampled in the data in order to be captured by the model. For example, in protein-ligand binding, a dissociation rate can only be computed if each step of the dissociation process has been sampled at least once.
Orders of magnitude of speedup can be achieved with enhanced sampling methods such as umbrella sampling, replica exchange dynamics, or meta-dynamics [20][21][22][23]. The speedup is achieved by coupling the unbiased ensemble of interest with ensembles at higher temperature at which the rare-events occur more frequently or by using biasing potentials allowing to "drag" the system across an energy barrier. With such approaches, accurate equilibrium expectations, such as free energy profiles can be computed efficiently, but the dynamical properties of the unbiased ensemble, such as transition rates, relaxation time-scales and transition pathways, are generally not available.
A common approach to reconstruct the kinetics from the free energy profiles is to employ rate theories such as transition state theory, Kramer's or Smoluchowski/Langevin models [24][25][26][27]. Such dynamical models introduce additional assumptions that cannot be selfconsistently validated because the predicted dynamics is not present in the data. Key assumptions underlying such rate models are usually the existence of a time-scale separation and approximate Markovianity on a single or few reaction coordinates -assumptions that are unlikely to hold for complex multi-state systems describing, e.g. bio-molecular dynamics.
Computation of kinetic quantities without rate models is possible with path sampling methods, such as transition path sampling [28], milestoning [29], transition interface sampling [30], and multi-state transition interface sampling [31]. A challenge is that in all of these approaches, the transition end-states must be defined apriori and all relevant rare-events must be distinguishable in the reaction coordinates, cores, or milestones that the method operates on.
Here we construct a general simulation approach that enables the computation of kinetic observables related to slow processes without having to explicitly sample the rare-events. It is based upon a very simple but general idea: Simulations are often constructed in such a way that they obey microscopic reversibility or at least a generalization thereof [32,33]. In this case, for any partition of state space into sets i, j, . . . and any choice of the lagtime τ , we have the detailed balance relation where p ij (τ ) is the probability of making a transition from set i to set j within a time τ and π i is the equilibrium probability of set i. Suppose we have knowledge about the equilibrium probabilities π i , π j from an enhanced sampling simulation. Then only the larger one of the two transition probabilities -p ij or p ji -needs to be sampled while the less probable event can be reconstructed by (1).
While this inference is trivial for a two-state system where three of the four components in (1) are known exactly, it is far from trivial for a system with many states and when some or all estimates are subject to statistical uncertainty. Here we establish a systematic inference scheme for combining multistate estimates of the equilibrium probabilities (π i ) with sampling data of at least the "down-hill" transition probabilities p ij . Our approach is built upon the framework of reversible Markov models [34,35] where (1) is enforced between all pairs of states. As a consequence, our estimates do not invoke any additional dynamical model, are accurate within a suitable state space discretization [8,11], and are precise in the limit of sufficient sampling.
The estimation procedure can reduce the sampling problem tremendously for processes with some long-lived states and some other states from which the system relaxes rapidly. This case is ubiquitous in meta-stable systems, because long-lived states are connected by shortlived transition states. But even long-lived states usually have very different lifetimes: For example many ligands or inhibitors bind to their protein receptor with nanomolar concentrations, meaning that the transition probabilities leading to the associated state are orders of magnitude higher than the dissociation probabilities. The present reversible Markov model approach lays the basis for estimating the kinetics and mechanisms of proteindrug dissociation by combining the much more rapid association trajectories with suitable enhanced sampling methods such as Hamiltonian replica-exchange [36] or umbrella sampling [18,37].

II. THEORY
A. Markov state models Classical dynamics, governed by Newton's equations in the case of an isolated system and by Langevin equation's for systems at constant temperature [38], gives rise to a transfer operator P propagating a phase-space density from time t to time t + ∆t [39][40][41]. Numerical solutions for Newton's or Langevin equations can be obtained for complex systems with many degrees of freedom, but a direct numerical assessment of the transfer operator is in most cases prohibited due to the curse of dimensionality.
Markov state models (MSMs) bridge this gap estimating the transfer operator on a suitably defined state space partition Ω = {s 1 , . . . , s n }.
(2) using trajectories obtained by direct numerical simulation [8,11]. MSMs model the jump process between states of this partition by a Markov chain. Observed transitions between pairs of states i and j are collected in a count matrix C = (c ij ) and the likelihood for the observed counts for a given transition matrix P = (p ij ) is given by While the likelihood functions allows to determine the maximum likelihood estimatorP optimizing the likelihood function for a given observation C over the set of all possible models P it does not specify the uncertainty of a chosen model.
For a finite amount of observation data there will in general be a whole ensemble of models compatible with the given data. In order to specify uncertainties and determine statistical errors of estimated quantities we need to infer the posterior probability of a model for a given observation. An application of Bayes' formula yields For a uniform prior, i.e. no a priori knowledge about the model, the posterior probability is given as a product of Dirichlet distributions B. Inference using a given stationary vector There are many methods that allow to efficiently estimate the stationary vector, even in situations in which a direct estimation from a finite observation of the Markov chain is infeasible due to the meta-stable nature of the system [20-23, 42, 43]. In such situations it is often possible to alter the system dynamics in a controlled way such that the artificial dynamics equilibrates more rapidly than the original one. The desired stationary vector of the original dynamics can then be related to the stationary vector estimated from the altered process [44][45][46][47][48].
In the following we want to show how such prior knowledge about the stationary vector can be used to improve the estimates of kinetic observables in systems with rareevents.
We are again given a finite observation of a Markov chain in terms of the count matrix C. Assume we are additionally given the stationary vector π for our system of interest and we know that the transition probabilities of the chain fulfill detailed balance for the given stationary vector, Then we can express the posterior probability for our model via (4). Prior knowledge about the stationary vector π in combination with the detailed balance assumption formally entails the following prior distribution on the posterior ensemble, According to (4) the constrained posterior is P(P |C, π) ∝ P(P |π)P(C|P ).
The effect of the prior (7) is a restriction of the posterior to the subspace of transition matrices fulfilling detailed balance with respect to the fixed stationary vector π.

C. Maximum likelihood estimate given a stationary vector
We can also use prior knowledge of the stationary vector to constrain the maximum likelihood estimateP to the set of matrices obeying (6) for a given stationary vector π. This results in the following convex constrained optimization problem which can be solved using methods outlined in [49].

D. Inference using a stationary vector with uncertainty
A stationary vector estimate usually carries a finite sampling error which should be accounted for when inferring a reversible transition matrix from data. From a Bayesian viewpoint we have to combine two sources of evidence. The observed count-matrix C from standard equilibrium simulations and the data from enhanced or biased sampling methods E used to estimate the stationary vector.
The posterior for transition matrices under the combined evidence P(P |C, E) can be decomposed as follows P(P |C, E) = dπ P(P |C, π, E)P(π|E, C).
The two approximations P(P |C, π, E) ≈ P(P |C, π) and P(π|E, C) ≈ P(π|E) lead to P(P |C, π) is the constrained posterior (8) and P(π|E) denotes an error-model for the stationary distribution estimated from enhanced or biased sampling simulations. Approximate sampling from P(P |C, E) can now be achieved by drawing a random sample π (1) , . . . , π (M) distributed according to a given error-model, π (k) ∼ P(π|E) and generating a sample of transition matrices P will then be approximately distributed according to P(P |C, E).
In [35] we have presented a Markov chain Monte Carlo approach to sample reversible transition matrices fulfilling detailed balance with respect to a fixed stationary vector. This method however has suffered from poor acceptance probabilities. In [49] we outline a method to efficiently generate samples from the constrained posterior using a Gibbs sampling algorithm that will be used here.
For given vector (π i ) detailed balance enforces a linear dependence between the transition matrix element p ij and the element p ji . As an immediate consequence the standard error of both elements for a sample generated from the posterior P(P |C) has to be equal, We will show how this can be used in order to significantly improve various estimates in situations in which p ij ≪ p ji .

III. RESULTS
In the following we will demonstrate the usefulness of (11) via a comparison of the standard error for kinetic quantities depending on rare-events that are either estimated from a Markov model of the direct unbiased simulation (unconstrained posterior (5)), as well as when enhanced sampling data is additionally used (constrained posterior (8) or constrained posterior with uncertain stationary vector (10)).

A. Finite state space Markov chain
Consider a three-state Markov chain with the following transition matrix The parameter b > 0 can be thought of as the height of an energy barrier between states one and three. The corresponding stationary distribution is given by The pair (π, P ) satisfies the detailed balance equation (6). Any process starting in state one has an exponential small probability of crossing over to state three. In fact a chain starting in state one can reach state three only via state two, but the probability to go from state one to state two is exponentially small in the barrier height b. The reversed process, going from state two to state one, occurs much faster. The same applies to state three and state two. The eigenvalues of this matrix are and the slowest time-scale in the system is given by, It is apparent from t 2 ≈ p −1 12 that estimates of t 2 and of p 12 have similar standard errors. The standard error ǫ for a matrix-element p ij for sampling from the unconstrained posterior (5) is For b = 4 and a single chain of length N ≈ 7 · 10 4 steps starting in state one we can on average expect c 12 = 4 resulting in a relative standard error of 50%. In order to decrease the error down to 1% we would need to run a chain of length N ≈ 100 2 ·10 4 = 10 8 steps. This is clearly an unsatisfactory situation and we would like to reduce the required simulation effort to reach a given error level as much as possible.
In comparison for an ensemble of M short chains of length L, with L ≪ 10 b , starting in state two one will on average observe a transition from state two to state one for every second chain, c 21 = M/2, so that a relative error of 1% for p 21 can already be achieved for M ≈ 10 4 , with L ≪ 10 b , so that the total simulation effort can be reduced by orders of magnitude.
We do not have explicit expressions for the standard error of matrix elements p ij when sampling from the restricted ensemble enforcing detailed balance with respect to a given stationary vector. It is however conceivable that the standard errors of p 21 can be reduced in the same way. The relation (11) guarantees that a small error for p 21 will also result in a small error for the rare-event quantity p 12 . Figure 1 shows the standard error of t 2 versus the total simulation effort. The error for a single long chain is estimated from a sample of transition matrices generated from the unconstrained posterior. The error for the ensemble of short chains is estimated from a sample of transition matrices generated from the constrained posterior using the algorithm outlined in [49]. From Figure 1 it is apparent that using a-priori information about the stationary vector in combination with an ensemble of short simulations started from the unstable state results in a three orders of magnitude smaller simulation effort when trying to estimate t 2 with a prescribed error. In particular, estimation of the rare-event kinetics can be conducted orders of magnitude before a direct simulation would even encounter a single transition event.

B. Doublewell potential
Let us now go to an example where the Markov state model is an approximation of the true dynamics. We employ Brownian dynamics in a double-well potential defined by The two minima of the potential at ±σ are separated by a maximum at −δσ/4, cf. Figure 2. The dynamics is governed by the following SDE, with dW t denoting the increments of the Wiener-process. The inverse temperature β = (k B T ) −1 controls the intensity of the stochastic fluctuations.
(15) defines a process where X t sample from the canonical distribution, The temperature dependent constant Z(β) is the partition function ensuring correct normalization, dx π(x) = 1. Spectral properties of this Markov process, such as the largest implied time-scale can be computed from a spatial discretisation of its associated transition kernel, cf. section A.   Figure 2. Potential V (x) and stationary distribution π(x),for Brownian dynamics in doublewell potential. The stationary distribution (shaded area) is scaled to fit the scale of the potential function. It can be seen that the stationary probability is concentrated in the meta-stable regions around the two minima of the potential at ±σ.
From an implied time scale estimation using a long trajectory with N = 5 · 10 6 steps we obtain a lag-time of τ = 10dt.
The stationary vector is either computed via numerical integration of the stationary distribution, (16) or estimated from umbrella sampling simulations using the umbrella integration method [50].
Estimates were computed using n u = 12 umbrella sampling simulations with M = 2 · 10 4 points per umbrella as well as from n u = 32 umbrella sampling simulations with M = 1 · 10 5 points per umbrella. To account for the uncertainty in the estimated stationary vector we used jackknife re-sampling [51] of the generated data and computed the stationary distribution for each re-sampled data set to model the ensemble of stationary vectors compatible with the observed umbrella sampling data.
In Figure 3 we show mean and standard error of the largest implied time scale t 2 versus the total simulation effort N . The total simulation effort N is composed of the simulation effort spent on obtaining a count matrix from standard simulations, N C , and the simulation effort spent on obtaining the stationary distribution from umbrella sampling simulations, N π , We compare four different approaches to estimating mean and standard error of the largest implied time-scale t 2 .
1. Generate a single trajectory starting in one of the meta-stable regions and compute estimates without a priori knowledge of the stationary vector.
2. Generate an ensemble of short trajectories starting on the barrier and compute estimates using the true stationary vector as prior information.
3. Generate an ensemble of short trajectories starting on the barrier and compute estimates with an error-model for the stationary vector as prior information.
4. Balanced sampling: split the total simulation effort equally between umbrella simulations and short trajectories starting on the barrier, N π = N C = N/2. Compute estimates updating the error model for the stationary vector according to the increasing amount of data available for the estimation.
Transition matrices are sampled according to (5) if no prior knowledge about the stationary vector is available, according to (8) if the true stationary vector is known a priori, and from (10) if the stationary vector was estimated from umbrella simulation data.
For the third approach we have estimated the stationary vector from a small as well as for a large amount of umbrella sampling data in order to demonstrate the dependence of the standard error of the kinetic observable on the error in the ensemble of input stationary distributions.
It can be seen from Figure 3 that for a fixed effort N π the standard error can not be reduced below a certain amount with increasing N C . This is a result of the nonzero statistical error in the estimate of the stationary vector for fixed N π . The usual N − 1 2 dependence of the standard error can be recovered for the proposed splitting N π = N C = N/2. Figure 3 shows the favorable scaling coefficient of such an approach leading to an almost two orders of magnitude faster convergence of the estimated quantity compared to using standard simulations alone.

C. Alanine-dipeptide
As an example for a rare-event quantity in a molecular system we use the mean first-passage time for the C 5 to C ax 7 transition in the alanine-dipeptide molecule. Alanine-dipeptide has been the long-serving laboratory rat of molecular dynamics [52][53][54][55][56]. The φ and ψ dihedral angles have been identified as the two relevant coordinates for the slowest kinetic processes of the system in equilibrium. The potential of mean force for the two dihedral angles is shown in Figure 4.
One can identify five meta-stable regions in the freeenergy landscape. The C 5 and P II regions correspond to dihedral angles found in a beta-sheet conformation, the α R and α L regions correspond to a right, respectively left-handed α-helix conformation. Reference values for the mean-first passage times between all pairs of sets have been computed from the maximum likelihood estimator of (3) using a total of 10µs of simulation data. Values can be found in Table I. For details of the computation of mean first-passage times see section B.
All computations were carried out on highperformance GPU cards using the OpenMM simulation package [57]. The used forcefield was amber99sb-ildn [58]  The latter approach allows to obtain a reliable estimate already before the average waiting time for a single rare-event τAB + τBA has elapsed. A comparison of the standard error b) indicates a two orders of magnitude speedup when estimating the rare-event sensitive quantity t2. The finite error for the estimate of the stationary vector results in a saturation of the error of t2 which can be further decreased using a more precise estimate of the stationary vector from additional enhanced sampling simulations. and the used water-model was tip3p [59]. The peptide was simulated in a cubic box of 2.7nm length including 652 solvent molecules. Langevin equations were integrated at T = 300K using a time-step dt of 2f s. The potential used for umbrella sampling simulations was In Figure 5 we show the estimate of the mean firstpassage time τ AB between the C 5 and the C ax 7 region together with the corresponding standard error ǫ(τ AB ) for different values of the total simulation effort N . The simulation setup is similar to the one described for the doublewell potential in the previous section. Instead of starting short trajectories directly on the barrier we start  Figure 4. Free energy profile of alanine-dipeptide as a function of the dihedral angles. Energies are given in kJ/mol. The average thermal energy kBT at 300K is 2.493kJ/mol. One can identify five meta-stable sets on the dihedral angle torus, here indicated by black lines. There are three low energy (high probability) sets C5, PII and αR with φ < 0 and two high energy (low probability) sets αR and C ax 7 with φ > 0. them from the meta-stable α L region. Figure 5 shows that combining umbrella sampling data and short trajectories relaxing from a meta-stable region with low probability (high free-energy) towards a meta-stable state with high probability (low free-energy) is able to achieve a standard error with an order of magnitude less simulation effort compared to an ensemble of long trajectories. The observed 8-fold speedup is in very good agreement with the expected speedup given by with τ AB = 60ns the mfpt for the slow "up-hill" transition from C 5 to C ax 7 , τ BA = 1.5ns the mfpt for the fast "down-hill" relaxation from C ax 7 to C 5 , and L = 5ns the length of individual short trajectories.
The present approach of estimating rare-event kinetics is much more powerful than traditional rate theories because quantities that can be estimated can be much more complex than only rates. As a reversible Markov model is estimated, full mechanisms, such as the ensemble of transition pathways from one state to another state can be computed. To illustrate this we compute the committor The mean first-passage time τAB of the C5 to C ax 7 transition is used as an observable for a rare-event process. a) Convergence of the mean value is shown for a small number of long chains starting in the C5 region (red), an ensemble of short chains starting in the αL region combined with umbrella sampling simulations. b) The standard error indicates a more than one order of magnitude speedup when estimating the kinetic characteristic of a rare-event τAB using short trajectories in combination with umbrella sampling simulations compared to using long trajectories and no additional information about the stationary vector.
probability function, cf. section C, from C 5 to α L using both estimates. It is seen that information about the stationary vector results in nearly the same committor function as one estimated using an order of magnitude larger simulation effort.
The solvated alanine-dipeptide molecule at room temperature, a system with more than 6000 degrees of freedom exhibiting meta-stable processes leading to average wall-clock times of more than 6 hours for the simulation of a single rare-event using specialized software adapted to high performance GPU hardware, demonstrates the usefulness of the proposed method for the study of rareevents in complex molecular systems. For this system reference values had to be obtained using a large amount of simulation data requiring a total of 50 days of GPU time. The equilibrium distribution of the system was estimated from umbrella sampling simulations using the umbrella-  Figure 6. Forward committor q + (x) for transition from C5 to αL region. a) shows a non-reversible reference estimate for N = 10µs of simulation data. Dark contour lines indicate the free energy profile. b) shows the difference between the reference estimate and a non-reversible estimate for N = 1µs of simulation data. There is a large error in the transition region due to insufficient sampling in the short simulation. c) shows the distance for an estimate using a combination of umbrella sampling and standard simulation data with N = Nπ + NC = 960ns. There is no significant error in the transition region, the small error close to the second saddle is probably due to insufficient sampling of this region by the reference simulation.
integration method, generating reliable estimates of the canonical ensemble restricted to a subset of relevant reaction coordinates already for simulation data obtainable in a single day of GPU time. For this system the proposed method is able to reduce the required simulation effort for the reliable estimation of rare-event sensitive observables by almost an order of magnitude.

IV. CONCLUSION
We have described a principle that allows for the first time to estimate rare-event kinetics efficiently without employing a restrictive dynamical model. Our approach is applicable when the kinetic properties of interest can be computed from a Markov state model discretization of the system.
The key idea of the presented approach is to use enhanced sampling methods to obtain reliable estimates of the equilibrium distribution in combination with direct simulations of the fast downhill processes. These data are combined rigorously in a reversible Markov model. Our approach can deliver estimates of kinetic properties, including rates, passage times, but also complex quantities such as committor functions and transition path ensembles while achieving orders of magnitude speedup compared to a direct simulation.
The present approach will be efficient whenever the rare-event occurs between low-probability and highprobability states. A very important example of this class is computational drug design, where the binding of the drug compound occurs relatively fast [18], while the unbinding may be many orders of magnitude slower. Yet the unbinding is crucial because it has been shown to be critical for drug efficacy [60].
While the applications in the present paper have used reversible Markov model estimates in such a way that the enhanced sampling simulation and the unbiased "downhill" simulations visit the same state space, the principle here can be generalized beyond this case. States visited only in one but not in the other simulation can be modeled by appropriate uninformative priors on the respective variables, e.g. uniform prior in the equilibrium distributions of states not visited in an umbrella sampling simulation. Moreover, the present inference principle can be exploited in an adaptive sampling framework [61,62] to optimally distribute the computational effort between enhanced sampling and unbiased molecular dynamics simulations.

ACKNOWLEDGMENTS
We are grateful to Feliks Nüske for stimulating discussions. This work was funded by the Deutsche Forschungsgemeinschaft (DFG) grant NO825/3-1 (BTS) and a European Research Council (ERC) starting grant pcCell (FN). The solution of (15) with initial position X 0 = x 0 on [0, T ] is usually carried out by choosing a regular discretization of the time interval 0 = t 0 < t 1 < · · · < t N = T.
with ∆t = t k − t k−1 for all k = 1, . . . , N . The evolution of the stochastic process is then approximated by the following time-stepping scheme with X 0 = x 0 and η being a N (0, ∆t) distributed random variable. The time-stepping scheme (A1) is known as Euler method or Euler-Maruyama method, [63]. For this simple time-stepping scheme the transition kernel of the resulting Markov chain is given by with x = X t and y = X t+∆t . p ∆t (x, y) is a Gaussian distribution with mean µ = x − ∇V (x)∆t and variance σ 2 = 2∆t/β.
The transition probability P ∆t (B|A) between two sets A, B can be computed from P ∆t (B|A) = A dx π(x) B dy p ∆t (x, y) A dx π(x) .

(A3)
Choosing a L such that p ∆t (x, y) is effectively zero outside of [−L, L] we pick a spatial discretization − L = x 0 < x 1 < . . . < x N = L (A4) with a regular spacing ∆x = x k − x k−1 for k = 1, . . . , N such that p ∆t (x, y) and π(x) are approximately constant on sub-intervals S i = (x k , x k+1 ]. In this case we have We can approximate the matrix elements p ij = P (S j |S i ) as p ij ≈ p(x i , x j )∆x. and compute spectral properties from the matrix (p ij ) using standard eigenvalue solvers. The covered material can be found in many introductory books to stochastic processes, cf. [64].
For a stochastic process (X t ) on a state space Ω the first hitting time T B of a set B ⊆ Ω is defined as The mean first passage time τ x,B to the set B starting in state x ∈ Ω is the following expectation value For a Markov chain on a finite state space Ω = {1, . . . , n} with transition matrix (p x,y ) the mean first-passage time can be computed from the following system of equations, Assuming that the chain has stationary vector (µ x ) we define the mean first-passage time τ A,B from set A to set B as the µ-weighted average of all mean first-passage times to B when starting in a state x ∈ A, Computing the mean first-passage time between two sets for a Markov chain on a finite state space with given transition matrix thus amounts to finding the stationary vector together with the solution of a linear system of equations -both of which can be achieved using standard numerical linear algebra libraries.

Appendix C: Committor functions
Committor functions have been introduced in the context of Transition Path Theory [12] and are a central object for the characterization of transition processes between two meta-stable sets.
Let (X t ) again be a stochastic process on a state space Ω and let A, B ⊆ Ω be two meta-stable sets. The forward committor q (+) (x) is the probability that the process starting in x will reach the set B first, rather than the set A, q (+) (x) = P x (T A < T B ). (C1) Again T S denotes the first hitting time of a set S. For a Markov chain on a finite state space with transition matrix P the forward committor solves the following boundary value problem [13], (C2)