Transition rate theory, spectral analysis, and reactive paths

The kinetics of a dynamical system dominated by two metastable states is examined from the perspective of the activated-dynamics reactive flux formalism, Markov state eigenvalue spectral decomposition, and committor-based transition path theory. Analysis shows that the different theoretical formulations are consistent, clarifying the significance of the inherent microscopic lag-times that are implicated, and that the most meaningful one-dimensional reaction coordinate in the region of the transition state is along the gradient of the committor in the multidimensional subspace of collective variables. It is shown that the familiar reactive flux activated dynamics formalism provides an effective route to calculate the transition rate in the case of a narrow sharp barrier but much less so in the case of a broad flat barrier. In this case, the standard reactive flux correlation function decays very slowly to the plateau value that corresponds to the transmission coefficient. Treating the committor function as a reaction coordinate does not alleviate all issues caused by the slow relaxation of the reactive flux correlation function. A more efficient activated dynamics simulation algorithm may be achieved from a modified reactive flux weighted by the committor. Simulation results on simple systems are used to illustrate the various conceptual points.


I. INTRODUCTION
The dynamics of a complex system with two long-lived metastable states is a classical problem that is often represented phenomenologically by the kinetic model, with the overall relaxation time τ * = (kAB + kBA) −1 and equilibrium probabilities p A = kBA/(kAB + kBA) and p B = kAB/(kAB + kBA).One of the most important theoretical frameworks to tackle such a problem in complex systems has been Chandler's activated-dynamics reactive flux formalism. 1Assuming that x is a good "reaction coordinate" for the system of interest, the forward transition rate kAB is then expressed as where k TST AB is the transition state theory (TST) rate evaluated at the transition state x † , [2][3][4] and the prefactor κ is the transmission coefficient serving as a correction (κ ≤ 1).In the context of this analysis, κ can be deduced from a time-correlation function evaluated at the molecular time scale, τm, expected to be much shorter than the overall relaxation time of the system τ * .For this reason, one of the most attractive features of the activated-dynamics algorithm is that κ can be determined from the fate of relatively short trajectories of length τm initiated at the transition state. 1,5,6enerally, to characterize the dynamics of complex systems with two long-lived metastable states, one may adopt two different fundamental perspectives.One may choose to focus on the global relaxation time τ * of the system or, alternatively, on the net unidirectional reactive flux J AB from the state A to the state B. For a simple two-state kinetic model, the global relaxation time τ * incorporates both the influence of forward A → B and backward B → A transitions, reflecting the natural back and forth dynamics occurring spontaneously in the system under equilibrium conditions.In contrast, the forward transition rate kAB, like the mean first passage time (MFPT), 7 is more directly associated with a specific unidirectional A → B reactive flux, picturing the kinetic of the system in terms of a "reactant" to a "product" state.In some sense, the first picture seems to simply reflect the system's natural kinetics, whereas the second picture frames the issue more deliberately by ascribing a directionality to the system's kinetics. 7Of course, both pictures are equivalent because under equilibrium conditions, the forward and backward net unidirectional reactive flux must be equal, J BA = J AB , and detailed balance implies that p A kAB = p B kBA.
Although this discussion is organized around an exceedingly simple two-state system, it helps illustrate the two different computational and theoretical framework that can be used to characterize the kinetics of highly complex systems.4][15] This construct, which aims its attention toward the set of A → B reactive trajectories forming the transition path ensemble, is a foundational element of the transition path sampling (TPS) algorithm [16][17][18][19][20] and transition path theory (TPT). 14While the spectral analysis is untainted by the somewhat subjective choice of the boundary states, these theoretical frameworks more directly focus on the specific transition that is the object of interest.][15] Originally introduced by Onsager in 1938, 21 and subsequently rediscovered in the late 1990s, [22][23][24] the committor probability is a fundamental building block of TPT. 14 This perspective is a critical insight in the formulation of the string method, 12,[25][26][27][28][29] aimed at determining the "reaction tube" that supports most of the reactive flux from A to B. It is the objective of the present analysis to shed new light on the characterization the kinetics of complex systems by examining the relationship between the reactive flux formalism, the spectral analysis of Markov models, and the TPT framework.

II. THEORETICAL DEVELOPMENTS A. Reactive flux formalism
For the sake of clarity, it is worth recalling the main elements of the reactive flux formalism for a system with two metastable states A and B. It is assumed that there is a single "reaction channel" between the two states.Having defined a population operator or indicator function HB that is equal to 1 when the system is in state B and equal to zero when the system is in state A, Chandler started by considering the normalized population time-correlation function, 1 where δHB(t) = HB(t) − ⟨HB⟩ is the deviation from the average.Additional useful relations include ⟨HB⟩ = p B , ⟨HA⟩ = ⟨1 − HB⟩ = p A , and ⟨δHBδHB⟩ = p A p B .Noting that C(t) ∼ e −t/τ * at long time according to the phenomenology, Chandler argued that the time derivative Ċ(t) must decay from an initial value to a plateau equal to −1/τ * in the limit of t → τm.In the context of this analysis, τm represents a molecular time scale that must be much shorter than the overall relaxation time of the system τ * = (kAB + kBA) −1 .Traditionally, one must then define a "reaction coordinate," x(z), as a function of a set of collective variables z.Once the transition state x † is identified, the indicator function HB can be constructed as where θ() is a Heaviside step-function.Assuming that x is a good reaction coordinate for the system of interest, the forward transition rate kAB = p B /τ * is then expressed as where x † is the position of the transition state, ρ eq (x) is the marginal equilibrium distribution along x, k TST AB is the transition state theory (TST) rate, and κ is the transmission coefficient serving as a correction to the transition state theory rate, Alternatively, the kAB transition rate can also be formulated as the time derivative of the conditional probability that the system will be found in state B at time t, assuming it was initially in state A at time t = 0, 30 = lim where ⟨ ḢB(t)⟩ = 0 was used.One of the most attractive features of the activated-dynamics algorithm based on Eq. ( 4) is that κ can be determined from the fate of relatively short trajectories of length τm initiated at the transition state x † . 1,5,6On the other hand, the approach may be The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpencountering considerable challenges when the time τm needed to reach the plateau of the time-correlation function of Eq. ( 5) is too long.While the formalism provides a sound description of the relaxation of a two-state system from a fundamental point of view, the relative measure of success or failure of the approach depends on the rate of convergence of Eq. ( 5).Methodologies that can help reduce τm and can lead to significantly better statistical properties of the transmission coefficient and the reactive flux estimates are important. 31,32

B. Effective propagator and spectral decomposition
The probability density of the system at time t is expressed as ρ(x, v; t), where x and v represent the set of coordinates xi and velocities vi, respectively. 12Using u ≡ (x, v) to represents a point in phase space, the forward propagation step (u → u ′ ) for the probability density from the time t to the time t + Δt is The elementary propagator for a null time step, P 0 (u ′ |u), is the identity δ(u ′ − u).The dynamical propagation, which we may formally represent as ρ(t + Δt) = P Δt ⋅ ρ(t), obeys the Chapman-Kolmogorov equation for arbitrary times.The propagator with the microscopic time step may be repeatedly applied an arbitrary number of times as ρ(t + nΔt) = (P Δt ) n ⋅ ρ(t).The forward-backward microscopic detailed balance relation,

Δt
is the backward propagator.While a formulation based on the microscopic propagator P Δt (u ′ |u) offers the most complete representation of the reactive paths, it is generally necessary to consider the dynamics projected onto a subspace of reduced dimensionality.We define the effective propagator Pτ for the finite lag-time τ within the subspace spanned by a subset of collective variables (CVs) as where z(x) = (z 1 (x), . . ., zN(x)) is a vector-valued function that maps every microscopic configuration x of the system on a set of values z(x).The reduced probability density of the system at time t is expressed as ρ(z; t).The forward propagation step (z → z ′ ) for the reduced probability density from the time t to the time t + τ is It is assumed that the dynamics within the reduced subspace of the CVs is Markovian with a finite lag-time τ and that the propagator obeys the Chapman-Kolmogorov equation, ρ(t + nτ) = (Pτ) n ⋅ ρ(t).It is assumed that the system is in equilibrium and that we have microscopic detailed balance, Pτ(z ′ |z) ρeq(z) = Pτ(z|z ′ ) ρeq(z ′ ), where ρeq(z) = ∫ dx dv δ(z(x) − z) ρeq(x, v) is the equilibrium probability in the subspace of the CVs.Under these conditions, the effective propagator Pτ(z ′ |z) yields a self-consistent representation of the dynamics of the system within this subspace (closure of the dynamical propagation), built on a trajectory generated via the elementary propagator P Δt (u ′ |u) with a time step Δt shorter than τ.
In practice, one should seek to determine the smallest possible lag-time τ that achieves Markovity for the effective propagator.An important framework to examine this issue is to rely on a spectral decomposition of the effective dynamical propagator. 8,9,12The righteigenvector ψ R k (z) of the operator is defined as 8,9 where the eigenvalue λ k (τ) = e −μ k τ .The constants μ k ≥ 0 represents the associated τ-independent intrinsic decay rate of the nth eigenmode.The eigenvector ψ R 1 (z) with the eigenvalue λ 1 = 1 (μ 0 = 0) corresponds to the invariant equilibrium vector, ρ eq (z).The eigenvalues are ordered from the slowest to the fastest process, i.e., 1 = λ 1 > λ 2 > λ 3 > ⋅ ⋅ ⋅, and 0 = μ 1 < μ 2 < μ 3 < ⋅ ⋅ ⋅.There is also a set of associated orthogonal left-eigenvectors, with and Orthonormalization can be expressed as The first right-eigenvector is actually the equilibrium distribution, ψ R 1 (z) = ρeq(z), and the first left-eigenvector is equal to unity.The equilibrium time-correlation function of an arbitrary function v(z) is where C. Spectral decomposition of a two state system Let us return to the normalized population time-correlation function between a "reactant" state A and a "product" state B of Eq. ( 2) that is the starting point of the activated-dynamics reactive The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpflux formalism. 1The indicator functions, HA and HB, are defined such that HA = 1 when the system is in state A and zero otherwise and that HB = 1 when the system is in state B and zero otherwise.By construction, HA + HB = 1, and the equilibrium probability of the state A and B is p A = ⟨HA⟩ and p B = ⟨HB⟩, respectively.Relying on a spectral analysis of the Markov dynamics, the long-time relaxation of the system displayed by the equilibrium time-correlation function of the indicator function HB for the state B, where If the indicator function HB matches the metastable states of the system accurately, the amplitude of the eigenvector ψ R 2 should be nearly constant for all the states i within the two metastable basins.The integral of the indicator functions with the other eigenvectors should be negligible, i.e., (HB ⋅ ψ R k ) ≈ 0 for k > 2. This means that the second eigenmode ψ R 2 essentially represents the global transfer of probability between the metastable basins A and B and the normalized population time-correlation function from Eq. ( 2) is Thus, if there is mainly one slow mode corresponding to the transitions between the metastable states A and B, the long-time relaxation of the correlation function is expected to reflect mainly the second eigenvalue μ 2 .In a system with two metastable states, μ 2 can be related to the overall relaxation time corresponding to τ * = 1/μ 2 and the forward and backward transition rates, kAB = p B μ 2 and kBA = p A μ 2 .If the indicator function overlaps with higher order eigenmodes, the dominant relaxation time can be revealed by considering the correlation function over a longer time τm = nτ, which is reasonable as long as e −μ 2 τ m ≈ 1 while e −μ 3 τ m ≪ 1.This is possible, in practice, only if there is a large gap between the eigenvalues μ 2 and μ 3 .

D. Perspective from transition path theory 1. Steady-state forward flux and committor probability
A spectral analysis may encounter issues when some of the slowest eigenmodes are uncorrelated with the A → B transition of interest.4][15] Within the TPT framework, 14 the net steadystate unidirectional reactive flux from A to B can be constructed from the probability to make a transition from z to z ′ multiplied by the probabilities that the trajectory arriving at z came from A and that the trajectory will then go on from z ′ to reach B, minus the probabilities that the trajectory arriving at z came from B and that the trajectory will then go on from z ′ to reach A. Invoking microscopic detailed balance, this yields where q(z), called the forward committor, is the probability that a trajectory started at z will first reach the state B. The bounds on the integral imply that the entire subspace of CVs has been divided in a region A ′ that includes the state A and a region B ′ that includes the state B. For the effective dynamics within the subspace z, q(z) must satisfy the backward propagation condition, with the constraints q = 0 when z ∈ A and q = 1 when z ∈ B. Equation ( 21) can be derived by counting only the forward flux z → z ′ arising exclusively from reactive A → B trajectories members of the transition path ensemble 14 or by considering the net unidirectional reactive flux from A to B under steady-state conditions. 13ecause one is free to move the boundary between the regions A ′ and B ′ , Eq. ( 21) can be transformed into a convenient unconstrained expression, 12,14,15 where we have defined the committor time-correlation function, Here, q(τ) is an implicit short-hand notation for q(z[x(τ)]).Equations ( 23) and ( 24) can also be used to variationally optimize The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpthe committor; 12,14,[33][34][35] minimizing J AB with respect to a trial function, i.e., δJ AB [q(z)]/δq(z) = 0, recovers the backward propagator condition of Eq. ( 22).Because the committor to reach state A before state B is simply equal to (1 − q), Eq. ( 23) makes it clear that the steady-state flux from state B to state A is equal to that from state A to state B, J BA = J AB .The forward transition rate can be determined as kAB = J AB /p A , and the mean first passage time (MFPT) from A to B is equal to 1/kAB.Additional features about reactive paths can be formulated.For example, the probability p AB (z) for a A → B reactive trajectory knowing that the system is at z is the probability 1 − q(z) that the trajectory came from A multiplied by the probability q(z) that the trajectory will lead to B. Averaging over all possible z weighted by the equilibrium distribution ρ eq (z) yields the total probability of a reactive A → B trajectory, p AB = ⟨q (1 − q)⟩, at equilibrium.The latter may be understood as the fraction of time a very long trajectory of length ttot spends on the transition pieces that are purely A → B reactive, τr/ttot.Since J AB is the total number of A → B reactive transitions nr during the time ttot, we can write ⟨q (1 − q)⟩ = (τr/nr) (nr/ttot) = ⟨τr⟩J AB .Thus, the mean transit time during a reactive transition ⟨τr⟩ is equal to ⟨q (1 − q)⟩/J AB . 13,20ecause J AB = J BA , then p AB = p BA and the mean transit time is independent of the direction.

Markovian dynamics in the short time limit
While the present analysis is strictly valid if the effective propagator within the reduced subspace of the CVs is Markovian with the finite lag-time τ, one may ask when the reactive flux J AB can be determined in the limit τ → 0, The flux is well-defined in the limit of τ → 0 if the underlying dynamics arises from a genuine continuous-time Markovian process, where K is the transition rate matrix, with Pτ ≡ e Kτ and P 0 (z Similarly, the flux is well-defined in the limit of τ → 0 if the effective dynamics within the subspace of CVs (z) is assumed to be diffusive with diffusion matrix D; it can be shown that (Appendix A) Under these conditions, the forward reactive flux J AB is given by the initial slope of the committor correlation function at τ = 0, i.e., the time derivative evaluated at τ = 0.If should be noted, however, that this expression would be identically equal to zero due to time reversibility if the microscopic dynamics were truly inertial. 36hysically, this suggest that the proper limit for τ must be larger than 0 to establish the diffusive dynamics, i.e., Ċqq(0 + ) ≈ (Cqq(τ) − Cqq(0))/τ, with a small but finite τ > 0.

Committor and eigenvectors
It is of interest to relate the unidirectional flux J AB to the overall relaxation time determined from the spectral analysis.For this, we need to have a model for the committor.While the committor q(z) is not quite a left-eigenvector of the effective propagator, a useful approximate construction can be written as 37 where The constants a and b must have opposite sign to guarantee that the vector is orthogonal to the equilibrium vector.It can be shown This construction makes a function q(z) that is equal to 0 for z ∈ A and equal to 1 for z ∈ B, which approximately satisfies the backward propagation condition of Eq. ( 22).Using the model committor, it can be shown (Appendix C) that the net unidirectional reactive flux from A to B is This expression is the steady-state flux from A to B for the simple two-state model.For this result to be valid, the lag-time τ must be sufficiently long to yield an effective propagation within the subspace z that is Markovian.If this condition is met, then the net forward flux J AB is independent of the value of the lag-time.If we assume that μ 2 τ ≪ 1, we have The expression shows that the steady-state flux from A to B remains the same if the lag-time is nτ, as long as nμ 2 τ ≪ 1.

Approximate committors
The formal analysis relies on the existence of the effective propagator Pτ(z ′ |z) representing the Markovian dynamics of the system within the subspace of the CVs, z.Nonetheless, it is understood that the microscopic trajectory x(t) is generated from the elementary propagator P Δt (u ′ |u) with a time step Δt that may be much shorter than the lag-time τ.Formally, one can always fall back on the elementary propagator to represent the time-evolution of the system.In this context, the broadest perspective is offered by treating the function q(z) as a function of the microscopic configuration x(t) that can be evaluated at any arbitrary time t, i.e., q(z[x(t)]).

ARTICLE scitation.org/journal/jcp
The time-correlation function Cqq(t) may be considered, even if the underlying dynamics is non-Markovian and comprises inertial aspects at times shorter than τ (this will be the case in Fig. 3).Obviously, the function q(z) associated with a Markovian dynamics within the subspace of the CVs is only an approximation to the committor q(x, v) in the full phase space.However, this construct does nonetheless provide a useful framework to analyze the true microscopic dynamics of the system.Importantly, the validity of the final rate expression does not require our representation of the effective dynamics within the subspace of the CVs to be exact at all times.A similar strategy was used by Ruiz-Montero et al. 32 to formulate an efficiency activated-dynamics reactive flux simulations for systems experiencing large dissipative forces (we return to this in Sec.IV).
If the dynamics of the CVs genuinely corresponds to a diffusion process or is governed by a continuous-time Markov process, then the committor time-correlation function Cqq(t) is rigorously linear for t > 0 + .In this case, there is no short transient at small times and the committor time-correlation function is immediately linear.The forward transition rate follows directly from the forward unidirectional flux J AB in Eq. (31).Whenever any of those conditions are not satisfied, there will be some transient behavior at short time and the variational forward flux relation will be valid only for a lagtime that is sufficient long to establish the linear behavior of Cqq(t).For example, the forward flux relation is not valid if the dynamical propagation within the subspace z is non-Markovian because the lag-time τ is too short or if an approximate committor function q ⋆ (z) is used.In this case, the equilibrium time-correlation function from Eq. ( 15) is where Such an approximate committor function may overlap with higher order eigenvectors, ψ L k with k > 2. The correlation function will relax to the dominant rate if some lag-time τm exists such that e −μ 3 τ m ≪ 1 while e −μ 2 τ m ≈ 1.In this case, the committor time-correlation function is not immediately linear, even if the dynamics within the subspace z is Markovian.As t increases, Cqq(t) should begin to vary linearly with time, a behavior that is valid only up to a certain time scale because Cqq(t) → ⟨(q⋆ − ⟨q⋆⟩) 2 ⟩ in the limit t → ∞.There is an intermediate finite time scale τm larger than the Markovian lag-time τ where the correlation function may be expected to vary linearly with time.Because of the expected linear dependence, the net forward unidirectional flux J AB is essentially given by the slope "s" of the correlation function Cqq(t).This argument is valid if the committor time-correlation function has the form Cqq(t) = c(t) + st, where c(t) is a small component that rapidly decays to some finite plateau value with lim Accordingly, we can write that the limiting slope of the committor time-correlation function is given by The microscopic time τm required for this relation to be valid may need to be longer than the lag-time τ that is needed to yield a Markovian dynamics within the subspace z.

III. ILLUSTRATIVE SIMULATIONS
For illustrative purposes, we consider three simple double-well one-dimensional systems with the potential, where z is in Å and W(z) is in kcal/mol.The potentials are shown in Fig. 1 (top).The symmetric wells yield equilibrium probabilities p A = p B = 0.5.In all cases, the double-well system is symmetric with a single free energy barrier, which, in turn, can be characterized as "narrow," "medium," and "broad."The Euler-Lagrange equation determining the committor probability q(z) is subject to the constraint q(z)| z≤z 1 = 0 and q(z)| z≥z 2 = 1.The solution is where k B T = 0.5915 kcal/mol.For the sake of simplicity, the boundaries z 1 and z 2 were set to −7 and +7, respectively, for the three cases.The calculated committors are shown in Fig. 1 (bottom).Each system was simulated with Brownian dynamics (BD), 30 assuming a diffusion coefficient D = 1 Å 2 ps −1 .A time step of 0.005 ps was used to generate a 1 μs trajectory (200 × 10 6 steps).Because the BD trajectory is Markovian along the z axis, the time step of 0.005 ps is the same as the lag-time τ.
The system with a broad barrier was also simulated with Langevin dynamics, 30 assuming a diffusion coefficient D = 1 Å 2 /ps and a mass m of 20 atomic-mass-unit.For these parameters, the relaxation time of the velocity-velocity correlation function k B T/Dm is 0.08 ps.At this time scale, the Langevin dynamics becomes essentially equivalent to a Brownian dynamics within a short time τ.A time step of 0.001 ps was used to generate a 1 μs trajectory (10 9 steps).
Three types of time-correlation functions were calculated from the BD trajectories.First, we consider the committor timecorrelation function Cqq(t) = ⟨q(0)q(0)⟩ − ⟨q(0)q(t)⟩, which was ARTICLE scitation.org/journal/jcpFIG. 1. Schematic illustration of a two-state system comprising a single narrow, medium, and broad free energy barrier.At the top (in blue), the potentials W(z) are defined via Eq.( 35) in kcal/mol as a function of z in Å.In the bottom (in red), the committors q(z) for the different potentials are calculated via Eq.( 37) with z 1 and z 2 set to −7 and +7, respectively.defined in Eq. (24).We also consider the indicator time-correlation function CBB(t) defined as The ratio in the parentheses is the normalized population timecorrelation function introduced in Eq. ( 2) to develop the reactive flux formulation.It is easy to show that the forward flux JAB = pA kAB = lim t→τ m ĊBB(t).Because CBB(t) is very similar in appearance to the committor-based time-correlation function Cqq(t), this form is more convenient to compare the two different formulations.Finally, we define the position time-correlation function as with δz(t) = z(t) − ⟨z⟩.By construction, we have CBB(t) = Czz(t) = Cqq(t) = 0 at t = 0, and all three time-correlation functions are expected to converge to a straight line with the same limiting slope as t increases (but remains much smaller than τ * ).The time derivative of Cqq(t) evaluated at the lag-time τ is related to the steady-state reactive flux, J AB according to Eq. ( 23).Following Eq. ( 4), the time derivative ĊBB(t) in the limit of t → τm also yields the steady-state reactive flux, although the molecular time scale τm may be longer than the lag-time τ.Finally, according to Eq. ( 15), the time derivative of Ċzz(t) is expected to behave similarly in the long time limit.

A. Reactive flux and spectral analysis
The activated-dynamics reactive flux time-correlation function of Eq. ( 4) must decay to a plateau value within a molecular time scale τm for the transition rate of the system to be well defined. 1,5,6Central to this argument is the assumption that the normalized population time-correlation function of Eq. ( 2) follows a simple exponential decay with the overall relaxation time τ * for times t ≫ τm.The spectral analysis of the effective propagator Pτ(z ′ |z) in terms of its eigenvalues μ k leading to Eq. ( 20) clarifies the context supporting this assumption.Relating the activated-dynamics to the spectral analysis by comparing Eqs. ( 2) and ( 20), we can see that to properly extract the overall relaxation time of the two-state system, τm must be chosen such that the conditions, μ 2 τm ≪ 1 and μ 3 τm ≫ 1, are satisfied.This is possible only if there exists a large gap in the eigenvalue spectrum, μ 2 ≪ μ 3 .This is a requirement that the two-state system must satisfy; otherwise, the theory based on the fluctuations expressed in Eq. ( 2) is inappropriate.
The spectral analysis also shed some light on important aspects of the activated-dynamics reactive flux algorithm to efficiently determine a transition rate from computer simulations.In practice, it is well understood that the algorithm is very efficient when the two states are separated by a sharp free energy barrier, the reactive flux time-correlation function in Eq. ( 4) rapidly within a relatively short time τm.However, it is also recognized that the algorithm may not be as effective in the case of slow diffusive motion over a broad flat free energy barrier.In this case, a much longer time τm may be needed to reach the plateau of the reactive flux timecorrelation function, requiring the simulation of long activated trajectories.

B. Illustrative simulations
To illustrate these points, Brownian dynamics simulations were carried out for simple one-dimensional systems (Fig. 1).The main results are shown in Fig. 2. The indicator time-correlation function CBB(t) is related to the activated-dynamics formalism to calculate the transition rate. 1 Position time-correlation functions, such as Czz(t), appear in the time-lagged independent component analysis (TICA). 38,39The committor time-correlation function Cqq(t) originate from the TPT framework. 7,12,14,15,33,40From Fig. 2, it is clear that all time-correlation functions ultimately converge to the same slope, reflecting the global relaxation time within the system.However, the different time-correlation functions reach the same slope differently.When the barrier is extremely narrow, the indicator time-correlation function CBB(t) converges extremely rapidly to a straight line with constant slope, whereas the position autocorrelation function Czz(t) takes a longer time.There are large fluctuations along z in the two broad wells, but those fluctuations are non-reactive-they are not associated with a A → B transition.In contrast, the fluctuations in CBB are faithfully reporting on such transition.In contrast, when the barrier is broad, the position autocorrelation function Czz(t) converges to a straight line more rapidly, whereas the indicator timecorrelation function CBB(t) takes a much a longer time.Nonetheless, even in the case of the broad barrier, the time-correlation functions relaxes to its constant slope within less than 1-2 ps.This is shorter than the time needed to diffuse from the center of the barrier to its edge, estimated as 2Dt = L 2 , with D = 1 Å 2 ps −1 and L = 7 Å, to be on the order of 25 ps.The latter time scale is comparable to the mean transit time during a reactive transition, ⟨τr⟩ = ⟨q (1 − q)⟩/J AB , 13,20 which is estimated to be on the order of 25 ps in this case.In the present simulations, J AB is estimated to be 4.4 × 10 −4 , 1.9 × 10 −4 , and 1.1 × 10 −4 ps −1 for the narrow, medium, and broad barrier, respectively.The corresponding global relaxation times τ * = (μ 2 ) −1 are estimated to be 568, 1315, and 2272 ps, which is much longer than the Markovian lag-time, τ (0.005 ps), and the mean transit time of reactive transitions, ⟨τr⟩.

C. Reactive flux and approximate committor
The activated-dynamics reactive flux formalism requires the definition of a dividing surface between the states A and B, whereas the TPT framework requires the definition of two boundary regions associated with the states A and B. Yet, it is interesting to note that the expressions derived above for the forward transition rate kAB from the activated-dynamics reactive flux formalism bears a certain similarity to that from the committor time-correlation function.Equations ( 4) and ( 34) are equivalent if one imagines that the indicator function HB is substituted for the committor q.Specifically, we have that Ċqq(τ) = ĊBB(τm) = JAB, with the substitution ⟨q(0)q(t)⟩ → ⟨HB(0) ḢB(t)⟩.Like the committor q, HB is equal to FIG. 2. Time-correlation functions calculated from Brownian dynamics simulations for the simple two-state system with a narrow, medium, and broad free energy barrier (shown in Fig. 1).Shown is the indicator autocorrelation function C BB (t) (in black), the committor correlation function Cqq(t) (in red), and the position autocorrelation function Czz(t) (in blue).The bottom plots show a zoom of the plots at the top.For a given potential, all the correlation functions converges to the same slope, which is related to the net steady-state unidirectional reactive flux J AB .For a narrow, medium, and broad free energy barrier, J AB is estimated to be 4.4 × 10 −4 , 1.9 × 10 −4 , and 1.1 × 10 −4 ps −1 , respectively.The tick marks along the y-axis are evenly separated by 0.005 in the top plot and by 0.0005 in the bottom plot (the y labels are not displayed for the sake of clarity).The time axis is in ps.

ARTICLE
scitation.org/journal/jcpzero in state A and one in state B, and in that sense, the Heaviside step-function HB stands as some kind of crude approximation to the exact committor.By extension, when there are multiple collective variables z, one could imagine defining the reaction coordinate as x = (z − z † ) ⋅ n and the indicator function as For HB(x) to serve as an effective (albeit crude) approximation to the committor q ⋆ (z), the vector n should be parallel to the gradient of the committor, ∇q(z † ), with z † in the vicinity of the transition state region.For example, z † can be defined as where q(z) = 0.5 is the separatrix.Because the indicator function in Eq. ( 40) serves as an approximate committor q ⋆ (z), a time τm longer than τ may be needed in this case to reach the linear regime of the time-correlation function.It is also interesting to consider Eq. ( 40) as a functional form depending parametrically on n and z † .Recalling that Eqs. ( 23) and ( 24) for the reactive flux J AB and time-correlation function Cqq(τ) provide a variational principle to determine the committor from a trial function, 12,14,[33][34][35] one would minimize the reactive flux J AB as a function of n and z † to obtain variationally optimized parameters.This sheds new light on the significance of variational TST, which seeks to optimize the reaction coordinate by minimizing the transition rate. 41It also suggests that the most meaningful reactive direction at the transition state in the multidimensional space of the CVs is parallel to ∇q(z † ), consistent with the multidimensional Kramers-Langer theory. 42Indeed, as noted in TPS studies, 17 configurations sampled according to δ (z − z † ) ⋅ n) ρeq(z) lead to a broad a distribution of committor probabilities when using a suboptimal reaction coordinate with a vector n that is not parallel to ∇q(z † ).
Approximating the committor as a Heaviside step-function is likely to be more accurate when a sharp free energy barrier separates the states A and B. However, in the case of diffusive dynamics occurring on the top of a wide flat barrier, such a Heaviside step-function is not a very good representation of the correct committor.This is precisely the situation where the activated-dynamics reactive flux algorithm encounters difficulties to converge rapidly (Fig. 2).This suggest that designing improved model committors may help formulate more efficient computational algorithms to determine the transition rate.
The idea of exploiting an approximate committor as a framework to improve the convergence of a transition rate calculation bares a certain similarity with a very interesting strategy previously designed by Ruiz-Montero, Frenkel, and Brey. 32Their method, which also starts from the idea that the Heaviside step functions in the high friction limit should be replaced with the smooth function akin to the committor, was shown to dramatically improve the statistical convergence of the reactive flux algorithm for diffusive barriers.While they used an initial perturbation and characteristic functions that resemble the steady state concentration profile, it is well known in TPT that the latter equal to ρeq(z)q(z)(1 − q(z)). 13 with the idea of using an approximate committor, the high friction steady-state picture served only to better formulate the problem in their method, and the validity of the final rate expression did not require this picture to be exact. 32

D. Committor as reaction coordinate
A striking observation from Fig. 2 is that the committor timecorrelation function Cqq(t) yields the correct slope at very short time in all cases, independent of the shape of the free energy barrier.This is in strong contrast with the time-correlation functions CBB(t) and Czz(t).The convergence of the committor time-correlation function to the correct slope within the lag-time τ is expected.If the underlying dynamics in the subspace z is governed by a continuous-time Markov process as in Eq. ( 28) or by a diffusion process as in Eq. (A1), then the committor time-correlation function is immediately linear and the transition rate can be determined in the limit of t → 0 + from Eq. (28).When those conditions are not quite satisfied, there will be some transient behavior over the time-lag τ needed to achieve Markovity for the effective propagator Pτ within the subspace z. 12 Establishing the Markovity of Pτ is generally non-trivial, although the present observations suggest that a simple criteria might be to verify whether the time derivative Ċqq(t) reaches a plateau as t → τ.Whether this conjecture can be demonstrated rigorously is unclear.Once this is done, the net steady-state reactive flux J AB and the transition rate kAB evaluated at the lag-time τ can be obtained on the basis of Eq. (23).
It is sometimes tempting to imagine that the dynamics along the committor is Markovian with the lag-time τ.To elaborate on this point, Cqq could be formally expressed as where Pτ(q ′ |q) represents an effective reduced propagator for the forward step (q → q ′ ) for the time step τ, which can be defined via Eqs.( 9) and (10).Let us recall that τ is the lag-time that achieves Markovity for the effective propagator Pτ within the subspace z. 12 However, there is no guaranty that the dynamics of the system projected onto the one-dimensional coordinate q ought to be Markovian with the same lag-time τ.Therefore, while Eq. ( 42) is rigorously correct, it does not imply that the effective propagator Pτ(q ′ |q) obeys the Chapman-Kolmogorov equation.][45][46][47] Nonetheless, there are reasons to believe that the committor might serve as a useful guide to define a reaction coordinate.Indeed, it has been shown that in the case of a multidimensional activated process controlled by diffusion, using the vector normal to the isocommittor plane separatrix (q † = 0.5) to construct a one-dimension reaction coordinate yields a rate constant that is identical to that predicted by the multidimensional Kramers-Langer theory. 42imilarly, as made clear by Eq. ( 40), the most effective 1D reaction coordinate in a multidimensional space z is parallel to the gradient of the committor in the region of the transition state, ∇q(z † ).However, simply treating the committor as a one-dimensional The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpreaction coordinate in the activated-dynamics formalism does not necessarily alleviate the issues caused by slow relaxation and intermediate time scales.In practice, redefining the population operator HB on the basis of the separatrix q † would return the same normalized population time-correlation function from Eq. ( 2) and the same reactive flux time-correlation function.In this case, the committor-based reactive flux of Eq. ( 45) will converge more rapidly than the activated-dynamics reactive flux time-correlation function Eq. ( 4).

E. Indicator-restricted committor time-correlation function
One of the most attractive features of the activated-dynamics reactive flux formalism as expressed by Eq. ( 4) is that it is based on a local sampling with trajectories initiated at the transition state, x † , along some one-dimensional reaction coordinate. 1,5,6This is in contrast with the reactive flux expression based on the committor time-correlation function Eq. ( 23), which is based on an unconstrained equilibrium average over all possible initial positions within the subspace z.Seeking to recover such a local sampling, we consider the steady-state reactive flux from A to B based on Eq. ( 21) expressed in terms of transitions from the point z with committor q(z) < q † to the point z ′ with committor q(z ′ ) > q † , This expression is formally correct for any value of 0 < q † < 1.
For the sake of convenience, we choose the separatrix where q † = 0.5.By identification, we can then write the state indicator functions as θ(q − q † ) = HA and θ(q † − q) = HB.We follow similar steps from Eq. ( 4), leading to the development of the activateddynamics reactive flux formalism.Defining the indicator-restricted committor time-correlation function, CAB-qq(t) = ⟨HA(0) HB(t) (q(t) − q(0))⟩, the reactive flux at the boundary between the A and B regions can be written as Like the activated-dynamics reactive flux formalism of Eq. ( 4), this expressions is based on a local sampling of trajectories initiated at a given interface, with q(z[x]) = q † .From this point of view, the expression focuses the calculation of the forward transition rate constant at a transition state.By analogy with activated-dynamics, one may note that the forward rate in Eq. ( 4) could also be expressed from the flux at the interface between the state A and B as kAB = ⟨HA(0) ḢB(τm)⟩/pA.However, it is noteworthy that this form is not equivalent to simply substituting x → q in Eq. ( 4).Equation ( 45) is not the same as simply using the committor as a "reaction coordinate" in the reactive flux formalism.What is different is the term [q(τ) − q(0)], weighting differently the dynamical excursion away from the dividing surface between the states A and B, which greatly affects the convergence of the reactive flux.This committor-based term is present here because the fundamental expression for the reactive flux at the boundary between the A and B regions Eq. ( 21) is derived by including only local transitions z → z ′ that belong to reactive paths with genuine A → B transitions. 12o illustrate these ideas and further examine the convergence of the activated-dynamics reactive flux formalism and transition state theory, 1,5,6 we consider again the time-correlation functions calculated from Langevin simulations in the case of the broad free energy barrier.The parameters were chosen such that the Langevin dynamics is in a high friction limit and becomes equivalent to Brownian dynamics within the short time-lag τ = 0.08 ps.The timecorrelation functions are shown in Fig. 3.The indicator function time-correlation function CBB(t) and its time-derivative ĊBB(t) corresponding to the reactive flux are shown at the top and bottom, respectively.The initial value of ĊBB(t) evaluated at t = 0 + corresponds to p A times the transition state rate k TST AB .The plateau value of ĊBB(τm) is equal to p A times the forward transition rate kAB, and the transmission coefficient κ(τm) is equal to ĊBB(τm)/ ĊBB(0 + ) via Eq.( 4).
The indicator-restricted committor time-correlation function CAB-qq(t) is shown in Fig. 3 (top).It essentially mirrors the evolution of the committor time-correlation function Cqq(t), rapidly converging toward the same slope.The reactive flux (bottom) shows that, like Ċqq(t), ĊAB-qq(t) reach a plateau within the lag-time τ of ∼0.08 ps (see the inset).In contrast, the activated-dynamics reactive flux ĊBB(t) decays to reach its plateau value after a considerably longer time τm (∼12 ps in Fig. 3).The position time-correlation function Ċzz(t) rises sharply and reaches its plateau faster than ĊBB(t) but still more slowly than Ċqq(t) and ĊAB-qq(t).The behavior of CAB-qq(t) displays some very interesting similarity with the TPS study of van Erp et al., 18 who showed that the time-derivative of the population correlation function constructed from indicator functions that depends on the history of the trajectories (past and future) also reached the correct plateau value at very short time.Whereas a trajectory's history (past and future) is explicitly visible in TPS, 16,17 such information in TPT is incorporated through the committor probability function. 14Accordingly, it is the term (q(τ) − q(0)) in the TPT-based expression for CAB-qq(t) of Eq. ( 45), which carries information about a trajectory's history and insures that only members of the transition paths ensemble are included, that dramatically affects the convergence of the reactive flux.While further analysis would be needed to formally relate the present analysis to TPS simulations, 18 this highlights the importance of attaching some information about the reactive or non-reactive "fate" of a given trajectory to determine a transition rate, a concept that has been well understood, from the early days of transition state theory (TST) [2][3][4] to more modern developments of the reactive flux formalism.The rapid convergence of the indicator-restricted committor time-correlation function CAB-qq(t) suggests that it may be possible to formulate an algorithm based on Eq. ( 45) that would rely on a local sampling like the activated-dynamics reactive flux formalism, 1,5,6 but which would converge more rapidly toward the correct transition rate constant.This advantage, however, needs to be qualified.Although such an algorithm would converge within a short lag-time, an accurate estimate of the committor is required.The latter can be calculated for a given configuration by shooting trajectories to determine the probability of reaching state B before state A, although the most practical approach is to optimize a trial function q(z) via the variational principle based on Eqs. ( 23) and (24). 12,34,35The important question is whether the usage of an approximate committor, which introduces additional factors associated high order eigenvalues as shown by Eq. ( 32), can be utilized productively.

V. CONCLUSION
The kinetic relaxation of a complex dynamical system dominated by two metastable states was examined from a Markov state eigenvalue spectral decomposition, [8][9][10][11] a committor-based steadystate reactive flux, 7,13,14 and the activated-dynamics reactive flux formalism. 1,5,6The present analysis was based on the assumption that the metastable states A and B were already known.9][50][51] Importantly, the robustness of the analysis and the overall shape of the committor function is typically not affected by the precise definition of the boundary states-this is one reason why the same definition was used for all three examples shown in Fig. 1 with very different free energy barriers.
The analysis showed that the different theoretical perspectives are consistent.This is most clearly displayed by the close correspondence between the eigenmode relaxation of Eq. ( 20), the steady-state reactive flux via Eq.( 31), and the activated-dynamics formalism via Eq.( 2).The spectral analysis of the effective propagator Pτ yields the eigenvalue λ k (τ) = e −μ k τ , with the second eigenvalue μ 2 corresponding to the global relaxation time τ * = (kAB + kBA) −1 = (μ 2 ) −1 of the two-state system, while the net steady-state reactive flux from A to B is J AB = p A kAB.
The significance of the different time scales that appear in the analysis of the kinetic relaxation time of a two-state system was clarified.While Δt is the microscopic time step used to generate the trajectory of the system, τ ≫ Δt is the shortest lag-time that achieves Markovity for the effective propagator Pτ within the subspace z.The analysis also revealed the existence of an intermediate time scale τm, which is associated with additional molecular processes that must decay to display the true global relaxation time of the two-state system.The intermediate time scale τm is larger than the lag-time τ but much smaller than the global relaxation time τ * .According to the spectral analysis, a good separation of time scales can be achieved if there is a large gap between the eigenvalues μ 2 and μ 3 such that e −μ 3 τ m ≪ 1 while e −μ 2 τ m ≈ 1.In particular, this intermediate time scale τm is displayed by the decay of the activated-dynamics reactive flux time-correlation function to a constant plateau corresponding to the transmission coefficient κ.However, it can also appear when an imperfect approximation for the committor is used in the expression for the steady-state reactive flux.In fact, there is a close correspondence between these two situations because the activated-dynamics reactive flux formalism is akin to using a Heaviside indicator function (population operator) to approximate the committor.This explains why the reactive flux with the indicator function is an effective route for a narrow sharp barrier but much less so in the case of a broad flat barrier.In this case, for example, the activated-dynamics reactive flux time-correlation function decays slowly to its plateau because the indicator function is affected by the diffusion process on top of the flat barrier.
An attractive feature of the activated-dynamics reactive flux formalism is the ability to calculate the forward transition rate constant from a local sampling of trajectories initiated at a given interface. 1,5,6While this provides an effective route to calculate the transition rate in the case of a narrow sharp barrier, it is much less effective in the case of a broad flat barrier.In contrast, the committor-based steady-state reactive flux rapidly converges to the The Journal of Chemical Physics ARTICLE scitation.org/journal/jcpcorrect value at the lag-time τ for both sharp and broad free energy barriers when the exact committor is used.The rapid convergence to the correct reactive flux within the time-lag τ and the absence of intermediate time scale τm is an attractive and robust feature of a committor-based steady-state reactive flux formulation.A tantalizing conjecture is whether this observation can serve to ascertain the Markovity of the effective propagator Pτ in the reduced subspace of the CVs.It was shown that the most effective choice to define a onedimensional reaction coordinate in the region of the transition state is along the gradient of the committor, ∇q(z), in the multidimensional subspace of collective variables z.However, while the committor can be exploited as a useful reaction coordinate, simply treating the committor as a one-dimensional reaction coordinate in the activated-dynamics formalism does not alleviate the issues caused by slow relaxation and intermediate time scales.In that sense, the idea that the committor is the "perfect" reaction coordinate, as is sometimes suggested, needs to be interpreted carefully.It is possible to express the reactive flux as a local average via the indicatorrestricted committor time-correlation function CAB-qq(t) in Eq. ( 45), which is one of the most attractive features of the activated-dynamics reactive flux formalism. 1,5,6Like Eq. ( 4), the forward transition rate constant can be calculated from a local sampling of very short trajectories initiated at a separating interface.
There is a large body of methods seeking to identify optimal reaction coordinates from complex dynamics, 11,33,47,52,53 including dynamical self-consistency, 46 memory reduction, 54 multi-dimensional spectral gap optimization of order parameters (SGOOP), 55,56 maximally predictive one-dimensional projection, 10 and variational committor-based steady-state reactive flux. 12,34,35Also related are the dynamical Galerkin approximation (DGA) formulated to predict the long-timescale behavior from short-trajectory 34,35 and nonparametric variational optimization of reaction coordinates. 57,58While challenging, the determination of an optimal reaction coordinate remains of great value, and a long term objective is to use these ideas to provide further guidance to theoretical frameworks built upon an optimized pathway. 29,59,60

FIG. 3 .
FIG.3.Time-correlation functions for the simple two-state system with a broad free energy barrier (shown in Fig.1).In the top, the indicator autocorrelation function C BB (t) (in black), the position autocorrelation function Czz(t) (in blue), the committor time-correlation function Cqq(t) (in red), and the indicator-restricted committor time-correlation function C AB-qq (t) (in green) is shown.The time axis is in ps.All the time-correlation functions converge to the same slope, which is related to the steady-state forward reactive flux J AB (top).The time-derivative of the time-correlation functions yielding the reactive flux are plotted (bottom) to illustrate how the slope relaxes to a constant plateau value corresponding to J AB = p A k AB = 1.05 × 10 −4 ps −1 , close to the value for the broad barrier obtained by Brownian dynamics reported in Fig.2.