PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptNIH Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Phys Chem B. Author manuscript; available in PMC Jul 26, 2013.
Published in final edited form as:
PMCID: PMC3406241
NIHMSID: NIHMS383241

Characterization of a Dynamic String Method for the Construction of Transition Pathways in Molecular Reactions

Abstract

We explore the theoretical foundation of different string methods used to find dominant reaction pathways in high-dimensional configuration spaces. Pathways are assessed by the amount of reactive flux they carry and by their orientation relative to the committor function. By examining the effects of transforming between different collective coordinates that span the same underlying space, we unmask artificial coordinate dependences in strings optimized to follow the free energy gradient. In contrast, strings optimized to follow the drift vector produce reaction pathways that are significantly less sensitive to reparameterizations of the collective coordinates. The differences in these paths arise because the drift vector depends on both the free energy gradient and the diffusion tensor of the coarse collective variables. Anisotropy and position dependence of diffusion tensors arise commonly in spaces of coarse variables, whose generally slow dynamics are obtained by nonlinear projections of the strongly coupled atomic motions. We show here that transition paths constructed to account for dynamics by following the drift vector will (to a close approximation) carry the maximum reactive flux both in systems with isotropic position dependent diffusion, and in systems with constant but anisotropic diffusion. We derive a simple method for calculating the committor function along paths that follow the reactive flux. Lastly, we provide guidance for the practical implementation of the dynamic string method.

Keywords: String method, transition path theory, drift vector, coordinate transformation, reaction mechanism

1. Introduction

The folding of proteins and their enzymatic catalysis are examples of reactions that proceed in a high-dimensional configuration space from an initial reactant state to a final product state. Reaction pathways aim to describe this evolution through configuration space and thus to capture the mechanism of the reaction. In molecular simulations, finding an ensemble of trajectories that travel along the entire reaction pathway from the reactant to product states is hampered by the long times needed to sample these generally rare transitions. One method for finding such transition trajectories is transition path sampling (see refs 1,2 and for more recent advances refs 35). Transition path sampling and more recently transition interface sampling6 help uncover the reactive trajectories connecting predefined reactant and product states through importance sampling techniques.1,2,6,7 These microscopic transition paths are expected to be highly diverse in full configuration space (e.g., with respect to the exact sequence of dihedral angle changes in a protein folding trajectory). However, one may hope that in the projection onto a suitably chosen set of coarse variables these transition paths are nearly coincident, and therefore a large fraction of these reactive trajectories are contained within a narrow tube around a dominant reaction pathway.

In a well-chosen space, such a pathway represents the mechanism of the reaction, in the sense that it describes the changes in key variables as the system progresses from a reactant to a product state along typical trajectories. One could identify the reaction pathway for this chosen set of coarse variables by projecting an ensemble of previously obtained transition paths from the full configuration space onto this coarse space. Alternatively, one can initially identify a set of coarse variables and then aim to directly (and more rapidly) sample the reaction pathway in the much smaller space of these coarse variables.

Coarse collective variables are typically chosen to capture the slow dynamics in (bio)molecular systems associated with reactive transitions. For example, typical solvent motions and bond vibrations are rapid compared to, e.g., torsion angles that depend on cooperative motion of multiple atoms. Global variables, such as those describing protein folding or large conformational changes, tend to relax even more slowly. Such coarse variables can be used as order parameters to report on the state of the system, and may also serve as reaction coordinates in simplified dynamic descriptions designed to capture the motions associated with large, slow, and rare transitions of macromolecules in a space of reduced and manageable dimensions.

The choice of coarse variables can have important consequences on both the convergence of pathway optimization algorithms and the interpretation of the results. Poorly chosen coordinates do not accurately parameterize the committor,2,8,9 which is defined as the probability for trajectories to proceed to a particular state, starting from a given point in configuration or phase space. As a result of the dependence on unresolved slowly relaxing degrees of freedom, estimators of both static and time dependent observables in the space of coarse variables will converge poorly. Among the consequences (and thus indicators!) of poor coordinate choices are hysteresis effects10,11 and large and often conflicting variations in the estimated local free energy gradients or drift vectors. Resulting estimates of free energy barriers tend to be inaccurate10,12,13 and misleading because the physically relevant barriers are not resolved. The assessment and optimization of reaction coordinates by exploiting the committor and related functions is thus the focus of intense ongoing efforts.9,1416,17,18,19

Here we explore different approaches for constructing dominant reaction pathways in a predefined space of coarse variables that are assumed to provide suitable reaction coordinates. These pathways provide one-dimensional representations of molecular reactions that serve as a basis for their mechanistic interpretation. In describing the slow dynamics in terms of transitions we assume that the configuration space can be divided into a set of identifiable states. For simplicity, we will focus on two states, reactant and product, but many of the concepts can be extended to systems with a larger number of identifiable states.

A variety of methods have been devised to find reaction pathways between two regions in configuration space. Searches in the space of paths aiming to maximize the flux20,21 or the path action22,23 result in boundary value equations where the reactant and product state must be defined a priori. The advantage of the action maximization is in resolving motion over arbitrary timescales, allowing one to identify motion also in the basins.22 However, the imposition of fixed end points constrains the evolution of the pathway and the results may be sensitive to the initial guess of the path.

As a possible alternative, the string method avoids such problems.24,8 Given molecular dynamics (MD) simulations in a large space (e.g., each atom’s coordinates) the string method provides a flexible framework for identifying configurations along the pathway from reactant to product states through a reduced space of coarse variables, and without the need to impose fixed end points. The pathway is represented by an ordered series of configurations connecting reactant and product regions. These configurations are iteratively updated to converge onto distinct final reaction pathways in the coarse variable space. One can use variants of the string method to find a minimum free energy path,25 a most probable path,26 a maximum flux path,27 or a principal curve (as in the finite temperature string (FTS) method28). As a major distinction, minimum free energy paths (MFEP)25 and principal curves28 follow the free energy surface and do not consider the influence of coarse variable dynamics on the reaction path, in contrast to the maximum flux path and the most probable path.

Our main focus is on pathway optimization methods that take into account the dynamics of sampling coarse variable space and aim to construct pathways that follow the reactive flux. The dynamics of the coarse variables is assumed to be effectively diffusive, arising from coupling between the underlying atomic coordinates, whether they obey Hamiltonian, Langevin or Brownian dynamics in the full space. Therefore, no assumptions need be made about the dynamics of the full coordinate space, as the string method is applied only to subsets of collective variables. An ensemble of multiple short MD trajectories, starting from a given initial condition in the space of coarse variables, can be used to determine the local drift vector of the coarse variable dynamics.29 These trajectories need to extend beyond a molecular time, in which the dynamics is dominated by fast inertial motions, to probe the relevant slow “diffusive” motion in coarse variable space.29 This estimate of the local drift has been used by Pan and Roux to update the string.26

There are two main reasons to use dynamic information in the construction of reaction pathways with string methods. First, the coarse-graining in time typically employed in the dynamic string method helps avoid getting trapped in local minima on a rough free energy surface in the coarse space. Second, local variations in the dynamics can result in favorable pathways redirected to pass through regions of higher diffusivity despite higher free energy. The idea of the minimum free energy path not necessarily carrying the highest number of reactive trajectories originated in early work by Berezkhovskii and Zitserman.30,31 These authors pioneered the concept of saddle-point avoidance, where reactive trajectories sacrifice free energy in exchange for faster mobility.

In many practical problems we expect substantial variations in the local diffusivity (or friction), even for a reduced set of Cartesian coordinates. One such example is protein folding, where the effective diffusivity of Cartesian-like coordinates such as the end-to-end distance was found to be one to two orders of magnitude smaller in the compact folded state than in the open unfolded state32. Relative motions of two groups on a protein chain tend to be dominated by solvent friction in the unfolded state, and by internal friction in the folded state, being hindered by the tight packing. In the case of protein folding, the assumption of position independent diffusion in the Cartesian sub-space of the end-to-end distance vector is thus clearly violated. Variations in diffusivity along coordinates also arise naturally in studies of stochastic dynamical systems.33

In this article we characterize the pathway constructed to incorporate dynamic information by following the direction of the average drift vector. We compare the flux and the committor q along this dynamic path (Figure 1), and distinguish it from the non-dynamic paths of FTS and MFEP methods. These comparisons are also relevant for related paths that were constructed without considering dynamic information, such as pathways determined using metadynamics approaches.34 To establish the properties of these paths even in the presence of non-constant, position dependent diffusion tensors we use coordinate transformations between different coarse-grained representations of the same dimensions (i.e., without averaging over coordinates). We start out in a system with constant diffusion, and perform a coordinate transformation to define a new coordinate system where variables will propagate with a position dependent diffusion tensor. Then, we compare the pathways sampled in each coordinate representation with the knowledge that the underlying physical systems are the same in both. This procedure helps unmask deficiencies in approaches in which the “optimum” path shows artificial dependences on the chosen coordinate system. We note that we are not comparing subspaces defined through different sets of coarse variables, but rather pathways obtained for the same subspace but with different coordinate representations related by linear or nonlinear transformations.

Figure 1
Schematic of dynamic string method. Multiple trajectories launched from the initial string images (blue circles) end at the black points. The averages of these end points determine the directions of the drift vectors (black arrows), and the new string ...

We first provide the theoretical background on the Smoluchowski equation for describing the dynamics of the coarse variables and on the committor function. Next we derive how the rules for transforming between two coordinate systems determine the relationship between the paths in the presence of general position dependent diffusion. We show that the drift path, unlike the MFEP or principal curve, follows the reactive flux within close approximation, and that the FTS path, which follows the gradient, depends on the coordinate representation of the system studied. We also derive an equation for calculating the committor along paths following the direction of the reactive flux. In the results section we provide numerical illustrations of the drift path for a selection of representative problems. Finally, we discuss the dynamic string method in more detail and the effect of parameter choice on the converged path.

2. Theory and Background

2.1. Diffusive Dynamics

2.1.1. Diffusion

We assume that a diffusion process accurately describes the dynamics in the space of coarse variables. Accordingly, the Smoluchowski equation defines the probability distribution p(x) of the N-dimensional coarse variable positions x = (x1,…, xN)T evolving in time

equation M1
(1)

where the diffusion tensor D(x) is a positive-definite symmetric matrix and in general a function of the coordinates, peq (x ) = Z−1 exp(−βV(x)) is the equilibrium distribution normalized by Z, V(x) is the potential, and β=(kBT)−1 is the inverse temperature T with kB Boltzmann’s constant. The derivatives are denoted by the column vector equation M2 and superscript T denotes the transpose. The dynamics of the coarse variables that sample this probability distribution are defined by the stochastic equations for (infinitesimal) time steps δt:

equation M3
(2)

The σ tensor35 can be obtained by a Cholesky decomposition, Dij = Σkσikσjk, and Rj are uncorrelated Gaussian random variable with zero mean and unit standard deviation. The last term defines the stochastic diffusive motion. The drift vector is defined by the first terms in parentheses,

equation M4
(3)

or in matrix notation, equation M5.

2.1.2. Committor

The committor function q(x) has a central role in the characterization of the reaction paths between reactant and product states.2,8 It defines the probability at each point that the system will reach the product state first before reaching the reactant state, and thus varies between zero and one. The surface of coordinates where q(x)= ½ (i.e., the isocommittor ½ surface) defines an ensemble of transition states. For the diffusive systems studied here, the committor36 satisfies the steady state backward Kolmogorov equation

equation M6
(4)

which is the adjoint of eq 1. The boundary conditions are q=0 in the reactant state, and q=1 in the product state. The quantity in parentheses is the drift vector component.

2.1.3. Reactive Flux

Another quantity central to the characterization of reactions is the reactive flux,37 equation M7 or in matrix notation,

equation M8
(5)

jR measures the flux of trajectories that reach the product state before returning to the reactant state, and as such is a subset8 of the total flux, defined as

equation M9
(6)

At equilibrium in a diffusive system described by eq 1 the total flux j is zero because of detailed balance, but the reactive flux jR is not as it travels in one direction. From eq 5 we see that the direction of the reactive flux vectors is determined by D (x)[partial differential]q/[partial differential]x. This implies that any path that follows the direction of the reactive flux also follows the direction of D(x)[partial differential]q/[partial differential]x. One implication is that for anisotropic D(x), the reactive flux will not be perpendicular to the isocommittor surfaces, as has been noted previously38 and is illustrated in Figure 2a. This means that the paths traveling along the reactive flux are not proceeding in the direction of the steepest increase in q to reach the product state. However, the speed of motion along coordinate directions is accounted for in the reactive flux direction.

Figure 2Figure 2
Illustrative examples of reactive fluxes, committors and reaction pathways. a) Directions of the reactive flux (black arrows) and normal to isocommittor surface (red arrow), with isocommittor contours from q ≈ 0 in brown to q ≈ 1 in tan, ...

2.2. Reaction Pathways

2.2.1. FTS and MFEP path

The FTS and MFEP are pathways obtained by the string method with distinct update rules. The FTS path is constructed under the assumption that in the low dimensional space of the chosen coarse variables, the dynamics, whether Langevin or Brownian, has a constant diffusion coefficient. One samples the equilibrium position of the coordinates in a constrained volume such as a Voronoi cell, and thereby follows the free energy gradient between the product and reactant well but with less localized sampling.28 The MFEP25 is constructed under the assumption that the dynamics in the Cartesian space of all positional coordinates was described by Langevin dynamics with a constant diffusion coefficient for all coordinates. Then, under the additional assumption that the chosen collective variables provide a good set of variables to describe the reaction, the collective variable dynamics should obey eq 2 with the diffusion tensor replaced by a metric tensor. This metric tensor accounts for the effects of the coordinate transformations on the displacement along each coordinate. The MFEP then follows the new free energy gradient, as transformed by the metric tensor. This method does not account for the dynamics of sampling the coarse variable space because the dynamics does not arise strictly from coordinate transformations, but depends in a generally complex nonlinear manner on the dynamics of the original coordinates.

2.2.2. Maximum Flux / Minimum Resistance Path

Berkowitz et al.20 showed that for a diffusive system with isotropic diffusion (Dij=D(x)δij), the path of maximum flux satisfies the coupled differential equations

equation M10
(7)

with s the arc length of the path and k the coordinate index. The derivation of this relation assumes that the minimum resistance path (MRP) x(s) in a Cartesian space follows the direction of the flux and that the flux is constant along the path. While the constant flux assumption is not generally true, it is nevertheless a reasonable approximation. Specific boundary conditions are also required in this derivation, such that the flux leaves the reactant state and arrives in a defined product state. The generalization of the MRP to systems with anisotropic and position dependent diffusion is shown in supporting information S3.

2.2.3. Drift Path

The drift path is defined by using the averaged drift vector of eq 3 in the coarse variable space. According to eq 2, for short times Δt an ensemble of trajectories starting at a point x(t) will on average drift by an amount Σj (−βDij[partial differential]V/[partial differential]xj + [partial differential]Dij/[partial differential]xjt in each direction i. The drift term in eq 2 can thus be estimated by monitoring the (deterministic!) displacement of the mean position29 from repeated simulations to average out the effect of random diffusive motion (or initial conditions in phase space conditioned on x). In contrast to the MFEP, the physical diffusion tensor entering the average drift describes the motion of the coarse variables as they sample their space. When diffusion is isotropic and independent of the position x, the drift vector will simply follow the direction of the gradient. For general position dependent diffusion, the FTS, MFEP and drift paths will differ. For details, see section 3.2.

For the case of isotropic diffusion, from eq 7 we see that the MRP follows the drift vectors if one ignores the curvature term D(x)d2 x/ds2. The curvature term does cause deviations between the MRP and the path following the drift vectors, because the MRP will not make turns as sharply as the drift vector. In Figure 2b we provide an example comparing the MRP (calculated using eq 7 above) with the drift path determined using the dynamic string method for a system with position dependent diffusion. The sources of deviation are confined to regions along the path where the drift vector is relatively small but sharply changing its orientation, such that the MRP will essentially cut the corners of the drift path to maintain a smoother path between the fixed reactant and product end points. This makes sense when considering that the drift vector depends only on the local energy surface and diffusion tensor, not on the fixed end points.

We therefore assume that the drift path largely follows the path of maximum flux for the isotropic diffusion case and below we will use coordinate transformations to extend this condition to anisotropic diffusion. Also, because the drift path is aligned with the flux, it will be roughly parallel to the reactive flux, assuming that distorting effects of the end-point constraints arising from the choice of reactant and product states can be neglected. Eq 5 thus also implies (assuming isotropic diffusion) that the drift path is roughly parallel to [partial differential]q(x)/[partial differential]x. This conclusion is consistent with previous work noting that in the subset of systems with constant diffusion [partial differential]V/[partial differential]x and [partial differential]q(x)/[partial differential]x are parallel,25,28 because for constant diffusion [partial differential]V/ [partial differential]x would also define the orientation of the drift vector. In Figure 2c we provide an example of the drift path for a rough28 Muller potential along with the isocommittor contours and the reactive flux vectors.

Thus far we have considered properties of the drift path for the case of isotropic and position independent diffusion. To examine how a non-constant diffusion tensor impacts the properties of reaction pathways, we use coordinate transformations between physically identical systems.

3. Effects of Dynamics on Reaction Pathways

3.1 Coordinate Transformation of the Smoluchowski Equation

With appropriately redefined potentials, corresponding equilibrium probability densities, and diffusion tensors, the Smoluchowski diffusion equation retains its functional form eq 1 in two different coordinate systems.39,40 When transforming from Cartesian coordinates xi to new coordinates zi(x1, x2, …, xN), the Jacobian matrix of the transformation, equation M11, and the matrix of the back transformation equation M12, account for volume changes under coordinate transforms, which here amount to entropic corrections. The fact that Jz (z) = Jx (x)−1 implies that equation M13 (with δik the Kronecker delta). The role of the Jacobian term in transformations to nonlinear coordinates has been discussed in several earlier works on free energy landscapes4143 and transition state theory, see e.g. 42,44,45. The supporting information section S1 details transformation rules for the different quantities entering the Smoluchowski equation. With eqs S2, S5 and S6, it is possible to determine how the committor, a path x(s) with s its arc length, the flux, and the drift vector will transform, both in orientation and in magnitude.

3.1.1 Flux Vector and Committor Transformations

Both the reactive flux and the full-time dependent flux change their orientation and magnitude under a coordinate transformation. Specifically, using the rules above we see that the flux transforms as

equation M14
(8)

The relative magnitude of the flux at different positions thus changes if the determinant |Jz (z)| of the Jacobian is not constant, which is the case for nonlinear coordinate transformations. Their orientation transforms as a tangent vector (discussed further below), leaving their direction invariant under coordinate changes. The committor function always transforms as a scalar (q(z) = q(x)x = x(z)).

3.2 Reaction Pathways Under Coordinate Transformations

Paths through coordinate space are one-dimensional curves that can be defined by the coordinates along the path and the tangent vector to the path. The path can be conveniently parameterized by its arc length along the curve s varying from zero to sf such that the tangent to the curve is a unit vector given by dx/ds. Under a transformation, the coordinates of the original path in x-space can be directly shifted to the corresponding coordinates in z-space. The path or curve is then redefined by transforming the tangent vectors at the corresponding coordinates in the z-space via dz = Jx (x)dx (i.e., these are contravariant vectors40). This transformation by Jx (x) is notable because it means that the vector in the z-space basis is equivalent to the original vector in the Cartesian basis.

Critically, because the original path transforms as a tangent, any path that follows the direction of the flux in x-space will continue to follow the flux in z-space, and vice versa. Another important implication of these transformation rules is that if a path follows the reactive flux in a system with arbitrary diffusion, it must follow the reactive flux and [partial differential]q/[partial differential]x in the back-transformed system with isotropic diffusion (assuming the back-transformation is possible40).

In the following, we study how the transformed original path relates to the new, z-space path constructed to follow the gradient or the drift vector in the z-space, and in which way the key properties of the new z-space path, where diffusion is not isotropic, differ from those of the original x-space path.

3.2.1 Gradient Path

If the original x-space path followed the gradient [partial differential]V/[partial differential]x, as in the FTS method, the corresponding transformed path will in general not align with the new gradient path [partial differential]Vz/[partial differential]z. This means that paths constructed to follow the derivative of the free energy will depend on the choice of coordinate system. The new z-space gradient is related to the original x-space gradient via

equation M15
(9)

Again for linear transformations the second term on the right-hand side (RHS) will be zero, but because the transformed gradient path is tangential to the vectors equation M16, the new path following the gradient according to eq 9 is clearly different. The gradient path transforms covariantly by the matrix Jz (z)T, in the same way as the basis vectors (eq S1), rather than by Jx (x), as a tangent vector would.

One implication of this difference in transformation rules is that along the new gradient descent path [partial differential]Vz/[partial differential]z is in general not parallel to [partial differential]q/[partial differential]z and therefore not perpendicular to the isocommittor surfaces. It would appear advantageous to construct paths that follow [partial differential]q/[partial differential]z to maximize the rate of increase in the probability of reaching the reactant state. However, in general the gradient path follows neither [partial differential]q/[partial differential]z nor the reactive flux for anisotropic diffusion tensors. Even under a simple rescaling of the coordinates, the gradient path is not invariant under coordinate transformations.

3.2.2 Drift Path

As an alternative to the gradient path, we now consider the drift path (Figure 1). The new drift vector in z-space is related to the original vector in x-space via

equation M17
(10)

This transformation rule applies if the original system has isotropic diffusion. In general, for a transformation from any arbitrary system with anisotropic position dependent diffusion, the components k of the second term on the RHS of eq 10 become40 equation M18.

For a linear transformation, Jx (x)is independent of position. As a result, the second term on the RHS of eq 10 vanishes. The transformed original path has new tangent vectors defined by Jx (x)drift(x)x = x(z). Therefore the new z-space drift path aligns exactly with the transformed original path. Because the flux transforms in the same way as the path (as a tangent vector), the new drift path will follow the flux just as the original path did in x-space. Also, because the flux changes magnitude uniformly everywhere (because |Jz (z)| is constant), the maximum flux path does not change and therefore the new drift path also continues to follow the maximum flux.

For nonlinear transformations, the new drift path is shifted from the transformed original path because of the second term on the RHS of eq 10. This means that for these transformations the new path is not guaranteed to be aligned with the reactive flux. However, in the examples considered below we find that the deviations tend to be small. Therefore, the transformed drift path should still follow closely the reactive flux. To justify this expectation, we note that the shift factor that produces the deviation (the second term on RHS of eq 10) is independent of the potential surface being sampled, whereas the drift vectors and the reactive flux are not. Hence only for shallow potentials where the path follows drift vectors of relatively small magnitude will the shift factor have a sizeable impact on the direction of the drift. In steeper potential valleys where the flux is more concentrated and the magnitude of the drift vectors is accordingly larger, this shift will have less effect on changing the drift path. For example, if we compared two systems with the same potential and diffusion constant but at different temperatures, at the lower temperature the drift vectors should be larger, with the impact of the shift therefore lessened, and the differences between the old and new paths should be smaller. At high temperature the paths will differ more, but the ensemble of reaction paths should then also widen, as the sampling of the potential becomes less sensitive to small variations in the surface. This is illustrated below in section 5.2.

Finally, under nonlinear transformations the Jacobian is not constant and the magnitude of the flux also changes. This means that although the transformed original path still follows the direction of the reactive flux, the length of the flux vectors have changed by |Jz (z)| and this could indicate that larger flux vectors are offset from the original path. Empirically, this change in magnitude seems to increase the size of the flux vectors in the direction opposite to the shift in the shift of the drift vector, suggesting that even if the new path aligns with the reactive flux, an optimized path would be displaced to follow the larger flux vectors. One hypothesis for an improved path is to follow the contravariant drift vector introduced by Graham39 instead of the normal drift. If the determinant of the diffusion tensor is independent of position, this amounts to following

equation M19
(11)

in z-space (otherwise, additional corrections enter). Eq 11 is further motivated by the observation that out of all the paths following the reactive flux (eq 5), assuming that variation in the magnitude of D [partial differential]q/[partial differential]z is small relative to variation in peq, the path with the highest probability defined by exp(−βVz(z)) should have the largest flux. Hence, this would be the path defined by D [partial differential]q/[partial differential]z that also follows the minimum of βVz(z), or the vector in eq 11. However, there is no simple algorithm for constructing this path, as multi-dimensional position dependent diffusion tensors are not trivial to extract from simulations, although examples4648 demonstrate it is certainly possible. Most directly, one can estimate D from the covariance matrix of the endpoints of the ensemble of trajectories run for a time Δt to construct the drift path, Dij ≈ cov(xi, xj)/2Δt.

4. Committor Along the Path

Along a converged drift path, we can estimate the value of the committor q using an analytical approximation. We first solve eq 5 for the gradient of q and project the gradient onto the reactive flux vector. An infinitesimal change of q along the path then becomes equation M20. Here we assumed that the drift path points in the direction of the reactive flux jR. If we further assume that the magnitude |jR| of the flux is constant along the path (as was assumed also in the construction of the MRP20), we obtain an approximate expression for the committor along the path,

equation M21
(12)

where D−1(s) = tTD−1(s)t with t a unit tangent vector of the drift path at arc-length s, and peq(s) [proportional, variant] exp(−βV(s)). Eq 12 generalizes earlier results,8 applies to any path that follows the reactive flux, such as the drift path, and is exact for one-dimensional systems (as has been previously noted9,49). The potential V(x(s)) = V(s) along the path can be obtained by integrating its gradient, i.e., the mean force, which can be estimated, e.g., from harmonically restrained simulations (see, e.g., ref 50). In supporting information section 2 we provide an alternative derivation of eq 12 starting from eq 4 and under slightly different assumptions.

In Figure 3 we compare the committors calculated along different paths using eq 12 with the full solutions to the 2D committor equation (eq 4). The agreement is very good, with only small deviations for the rough Muller potential, suggesting that the approximations made in the derivation of eq 12 are reasonable.

Figure 3
Committor along the path as a function of arc length s. a) For the potential shown in Figure 2a, the committor along the drift path is calculated using eq 12 (blue curve), and the solution to the full committor equation (eq 4) is shown at each point along ...

5. Results and Discussion

5.1. Linear Coordinate Transformations

We first study the double-well potential defined in ref 38 and pictured in Figure 4a. For this system, the MRP and the drift path are in excellent agreement. We then perform a linear transformation to a stretched coordinate system, z1=αx1 and z2=x2, with α=5 (Figure 4b). In this stretched coordinate system, the diffusion tensor is still diagonal but anisotropic, such that D11= α2, and D22=1. The determinant of the Jacobian is a scalar, α−1, and drops out of both the committor equation and the Smoluchowski equation. Hence the potential V sampled in z-space is not modified from the one sampled in x-space (except for a constant shift factor). As we noted above for linear transformations, the new drift path follows the transformed original path, and the magnitude of the flux changes uniformly everywhere. Therefore, the new drift path is parallel to the reactive flux direction D [partial differential]q/[partial differential]x and follows the maximum flux path (as it did in the original x-space).

Figure 4
2D double-well potential with isotropic and anisotropic diffusion. a) Isotropic diffusion with D=1. For this potential38 the drift path (blue) and the MRP (magenta) are in nearly exact agreement. The isocommittor surfaces capture the progress of the reaction ...

For completeness, we note that the MRP defined by the differential equation of Berkowitz et al.20 does not simply transform as a tangent vector, and therefore does not align exactly with the transformed original MRP. Since the flux itself does transform as a tangent and in this example |Jz (z)| is constant, this means the MRP is no longer in complete agreement with the maximum flux path. This inconsistency in the behavior of the MRP under transformation appears to be due to the curvature constraints on the path resulting from fixed boundary points and the assumption that the flux along the path is constant.

5.2. Nonlinear Coordinate Transformations

To demonstrate the properties of the drift path in a system with general position dependent diffusion, we start from an x-space system with constant diffusion sampling the Muller potential, shown in Figure 5a (defined, e.g., in ref 27). In one nonlinear transformation of the Muller potential, we define z1=(x+3)(y+1) and z2=y. This transformation results in a position dependent Jacobian and diffusion tensor in z-space. However, in this example (Figure 5b) the new drift path does in fact align exactly with the transformed original path because the second term on the RHS of eq 10 vanishes, equation M22. Therefore in this case the new drift path follows the reactive flux just as the original path did. However, the magnitude of the flux did change because |Jz (z)| is position dependent. So again, the larger flux seems to be shifted slightly away from the transformed drift path.

Figure 5Figure 5
Nonlinear transformations of the Muller potential. a) In the original coordinate system the simulations are run with diffusion constant D=1 and kBT=9. The MRP in magenta and the drift path in blue are in close agreement. b) After the nonlinear transformation ...

In another example, we perform the nonlinear transformation, z1=(x+2)2 and z2=(y+1)2 (with x>−2 and y>−1). In this case, from eq 10 we see that the direction of the new z-space drift path is shifted from the transformed original path by the constant vector equation M23 in the z1 and z2 directions. As we noted above, this shift factor is independent of the temperature. When we compare the new drift path and the transformed original path at a relatively low temperature kBT =9 creating a barrier of ~9kBT, they are nearly the same, and both closely follow the features of the potential surface (Figure 5c).

By scaling this Muller potential by a factor of 0.18, we create a much shallower potential to start with, such that the barrier is only ~1.7kBT (Figure 6a). For this system the direction of the reactive flux is less sensitive to the shallower potential surface and the MRP will actually pass through the saddle region at very slightly higher energy (even with constant diffusion) to maintain lower curvature of the path. After performing the same coordinate transformation as above, the new drift path is slightly offset compared to the transformed original path (Figure 6b). With the increased width of the high flux region, both paths appear to follow the reactive flux fairly accurately. However, regions of large reactive flux are shifted in the opposite direction compared to the new drift path. We note that this higher flux path appears to agree well with the path postulated above in eq 11.

Figure 6Figure 6
Shallow Muller potential. a) Results for a Muller potential scaled by a factor of 0.18 with a barrier of ~1.7kBT, D=1 and kBT=9. Because of the shallowness of the potential, the MRP (solid magenta), even in this constant diffusion case, cuts across the ...

5.3. Isotropic, Position Dependent Diffusion

In Figure 2b we used an example system where the diffusion tensor is position dependent but isotropic, as defined in ref 51. This system provides an example of saddle-point avoidance, where the paths carrying the maximum flux avoid the lowest energy region of the transition surface because the diffusion is slowed there. This diffusion tensor is imposed, rather than resulting from a coordinate transformation, such that a position independent diffusion cannot be obtained through a variable transformation.40 However, because the diffusion is isotropic, the drift path will still follow the maximum flux. The small deviations between the MRP and the drift path in Figure 2b are due to the lack of curvature constraints imposed on the drift path.

6. Implementation of Dynamic String Method

A possible implementation of the dynamic string method has been described elsewhere.26 Here we comment on the role of the algorithm parameters in defining the final path. To briefly summarize the method (see also Figure 1), an initial string is defined between proposed reactant and product states by connecting a series of configurations through, e.g., interpolation. Following ref 29 each point along the string is initialized by drawing all degrees of freedom from equilibrium distributions conditioned on the coarse variables (see Section 6.3). From such configurations, multiple trajectories are launched without constraints and run for a time Δt. The difference between the average final position and the initial position, divided by Δt, is used as an estimate of the drift vector. The points along the string are moved along these vectors to their new positions. The string is then reparameterized to maintain equal arc length segments between the points (or images) along the string. The string endpoints can be updated to refine their positions.

6.1. Convergence Parameters

The major parameters to choose when implementing the algorithm are the number of images along the string, the number of trajectories to run from each image, the number of iterations to converge to the final string, the time length Δt of the trajectories, and any smoothing or curvature constraints on the string shape. The first three parameters listed largely control the convergence of the string. Beyond a certain limit, additional images, trajectories, or iterations will not significantly alter the location or shape of the final string. If too few images along the string are used, the final path will capture only the major features defining the flux. The resolution can be improved by adding more images to a converged string. Too few trajectories will introduce noise into the drift vector estimates and slow down the convergence of the string.

6.2. Time Resolution of the Path

The time length Δt of the short MD trajectories is the most important parameter in defining the final drift string. Changing the time length of the drift trajectories can change the location of the final converged string. However, the effect of trajectory length on the final string is predictable. This dependence on the length of the drift is consistent with defining the mechanism of the reaction over a specific coarse-grained time resolution. Short trajectories resolve faster motions and therefore more features of the potential surface. For very long trajectories, the images along the string will have sufficient time to reach, or nearly reach, the minima on the free energy surface. As a result, all images will pull the string towards either the reactant or product state, and the reparameterized string will stretch linearly from one state to the other. For intermediate lengths this pulling will be less pronounced and more localized, allowing the string to cut corners in shallow regions of the energy surface. This smoothing or pulling of the string is similar to adding a curvature constraint on the strings’ shape.

In the numerical example systems, the true dynamics of the systems was diffusive, and as a result, the shortest possible trajectories still followed the direction of the drift vector. However, when running simulations from full molecular systems, trajectories that are too short will not allow the atomic degrees of freedom coupled to the collective variables enough time to relax, preventing the coarse variables from ‘drifting’ to their new position. Hence the trajectory lengths for biomolecular systems cannot be taken arbitrarily small, and Δt should be chosen appropriate for the system and coordinates to sample the relevant diffusive dynamics, not the underlying inertial dynamics.

6.3. Numerical Simulation Methods

Simulations of the example systems were run by propagating eq 3 with a time step of δt=10−4. Convergence of the string was checked by comparing the final string every 1000 iterations. Generally between 2000–5000 iterations were needed, depending on the potential (which were performed here without attempts to optimize the computational efficiency). The initial string was created as a line that did not need to start and end in the reactant and product state, as the algorithm will update these positions along the string. The drift vectors were collected from trajectories that ran for Δt =5–10δt from each image along the string. The number of trajectories launched from each image was ~20–30, although higher numbers were tried and results were not sensitive to the number of trajectories used. Around 30–50 images were used to determine the string, with more images being used for strings with longer arc length. In the particular examples, we found that including substantially more string images only introduced more noise in the sampling of the string rather than resolving any additional features of the energy surface. At every iteration of the string update, each image is moved along the direction of the drift vector a distance Δτ towards its new location. This stepsize Δτ is similar to an overrelaxation parameter and can be much larger than δt. It was set at values of ~0.05, although the results were not sensitive to this stepsize. Simulations with Δτ=0.01, or 0.1 produced the same path, although for smaller Δτ more iterations were needed.

The committor function (eq 4) was calculated numerically using the PDE toolbox of MATLAB. To define the reactant and product state for the committor calculation, small ovals were centered at the known basin minima. Neumann boundary conditions set to zero were used for the edges of the coordinate domain. The MRP (eqs 7 and S14) was determined using the boundary value solver in MATLAB. The reactive flux was calculated from the numerical derivative of the committor function using eq 5.

6.4. Initializing the String Images

In practice, the space spanned by x is only a subspace of the full configuration space, such that initial conditions have to be obtained by equilibration conditional on x, and trajectories have to be projected onto x. Before multiple trajectories can be launched from each image along the string, multiple configurations representing each image must be generated. This can be accomplished by running simulations at each image where the coarse variables are constrained (or restrained) to their current values, but the other degrees of freedom are sampled from a conditioned equilibrium distribution.

7. Concluding Remarks

Reaction pathways through the conformational space of complex (bio)molecular systems provide mechanistic insight into transitions of configurations and the rates of such transitions. The reactive flux provides a basis for defining a reasonable transition path that captures where most reactive trajectories will pass through configuration space. To construct pathways that follow the reactive flux, we advocate the use of the dynamic string method. For the cases of position independent anisotropic and position dependent isotropic diffusion the resulting paths follow the direction of the maximum reactive flux and are independent of the particular choice of coordinate representation in a given coarse variable space. Even for cases with general diffusion tensors, we find that the resulting string remains close to the direction of the reactive flux, with larger deviations expected for shallow potentials.

In contrast, reaction paths constructed to follow the free energy gradient vary artificially with the coordinate representation and follow the reactive flux only for systems with constant and isotropic diffusion. Paths constructed to follow the MRP would be an improvement over these free energy paths, but finding the MRP is problematic as it solves a boundary value problem that requires fixed end points and is sensitive to the roughness of the potential surface.

The practical implementation of the dynamic string method is straightforward. By following the averaged drift vector to update the string to convergence, the dynamic string method is able to avoid local small traps on the free energy surface, similar to the FTS method and in contrast to pure gradient based methods like the original string method and differential equation based methods21 like the MRP. It also does not require fixed reactant and product states but can move the string to locate these end points of the path. For relatively flat or shallow regions of the energy surface, the FTS method of sampling the equilibrium position results in noisy estimates of the local string position. The dynamic string method, however, reduces this noise by following the drift vector. By averaging the drift over many trajectories, small traps on the energy surface can be avoided. As a result, the path is more sensitive in locating end states within wide basins, but not so sensitive as to become trapped in local minima. Importantly, because the drift vector combines the effects of both the free energy gradient and the local diffusion tensor, the diffusion tensor does not need to be explicitly constructed from the simulations; in contrast, the MFEP requires the explicit construction of the metric tensor for the path updates. The dynamic string method does require choosing the length of the trajectories that define the drift vector and this can affect the shape of the converged path, as discussed in Section 6.2. However, this parameter can be interpreted as controlling the time resolution of the path in terms of the types of motion it captures, and is somewhat analogous to the choice of force constants needed to constrain images in the MFEP method.

The drift path is shaped by local features of the potential surface and unlike the reactive flux does not rely on defining specific reactant and product states. While this is an advantage in locating the path from unknown starting and end states, it excludes any curvature effects on the shape of the path. The effect of curvature in shaping the path is more apparent in regions where the path changes direction, and for potentials that are relatively shallow. A curvature constraint can be included in the updates of the string path, as described in ref 28, which will also help smooth the path.

We demonstrated that for any path following the reactive flux the committor along the path can be obtained by integration, eq 12. The behavior of the committor is dominated by the variation in the unscaled potential along the path, which appears in the exponential rather than as a scale factor like the diffusion tensor. The results for calculating the committor along the path do not apply to the minimum free energy paths or the FTS path, except in the case of constant diffusion, because in general these paths do not follow the reactive flux. Finally, we note that even paths constructed to follow the reactive flux do not necessarily provide an optimal one-dimensional reaction coordinate for calculating free energy barriers and transition rates when the reaction rate is controlled by travel through the saddle region. The gradient of this optimal one-dimensional reaction coordinate follows that of the committor, in particular near the transition point.38 In the case of anisotropic diffusion this direction is not coincident with the reactive flux or the free energy gradient, but can be estimated via eq 5 from the reactive flux and the local diffusion tensor.

In conclusion, pathways following the gradient of the free energy only account for structural or energetic changes in conformations, and are not invariant under simple reparameterizations of the collective variables used in their construction. As a practical alternative, the drift path is less sensitive to reparameterizations and, in addition, includes relevant dynamic information. Even if it does not perfectly track the reactive flux in all systems, the drift path captures the effects of both structure and mobility in the dominant reaction pathway.

Supplementary Material

1_si_001

Acknowledgments

We thank Dr. Attila Szabo for many stimulating discussions. This work was supported by the Intramural Research Program of the NIDDK, NIH. The research used the Biowulf Linux cluster at the NIH.

Footnotes

Supporting Information Available: The online supporting information contains three parts. A section (S1) defining how the terms entering the Smoluchowski equation change under a coordinate transformation, a section (S2) showing an alternate derivation of eq 12, and a section (S3) defining the MRP for systems with arbitrary diffusion. This material is available free of charge via the Internet at http://pubs.acs.org.

References

1. Dellago C, Bolhuis PG, Csajka FS, Chandler D. J Chem Phys. 1998;108:1964–1977.
2. Bolhuis PG, Chandler D, Dellago C, Geissler PL. Annu Rev Phys Chem. 2002;53:291–318. [PubMed]
3. Escobedo FA, Borrero EE, Araque JC. J Phys Cond Mat. 2009;21:333101. [PubMed]
4. Dellago C, Bolhuis PG. Adv Polym Sci. 2009;221:167–233.
5. Peters B. Mol Sim. 2010;36:1265–1281.
6. van Erp TS, Moroni D, Bolhuis PG. J Chem Phys. 2003;118:7762.
7. Pratt LR. J Chem Phys. 1986;85:5045–5048.
8. Vanden-Eijnden EEW. Annu Rev Phys Chem. 2010;61:391–420. [PubMed]
9. Best RB, Hummer G. Proc Nat Acad Sci USA. 2005;102:6732–6737. [PubMed]
10. van Erp TS. J Chem Phys. 2006;125:174106. [PubMed]
11. Rosta E, Woodcock HL, Brooks BR, Hummer G. J Comput Chem. 2009;30:1634–1641. [PMC free article] [PubMed]
12. Truhlar DG, Garrett BC. Annu Rev Phys Chem. 1984;35:159–189.
13. Bolhuis PG, Dellago C, Chandler D. Proc Nat Acad Sci USA. 2000;97:5877–5882. [PubMed]
14. EW, Ren W, Vanden-Eijnden E. Chem Phys Lett. 2005;413:242–247.
15. Peters B, Trout BL. J Chem Phys. 2006;125:054108. [PubMed]
16. Ma A, Dinner AR. J Phys Chem B. 2005;109:6769–6779. [PubMed]
17. Peters B. J Chem Phys. 2006;125:241101. [PubMed]
18. Chodera JD, Pande VS. Phys Rev Lett. 2011;107:098102. [PubMed]
19. Kirmizialtin S, Elber R. J Phys Chem A. 2011;115:6137–6148. [PMC free article] [PubMed]
20. Berkowitz M, Morgan JD, McCammon JA, Northrup SH. J Chem Phys. 1983;79:5563–5565.
21. Zhao R, Shen J, Skeel RD. J Chem Theor Comp. 2010;6:2411–2423. [PMC free article] [PubMed]
22. Olender R, Elber R. J Chem Phys. 1996;105:9299–9315.
23. Faccioli P, Sega M, Pederiva F, Orland H. Phys Rev Lett. 2006;97:108101. [PubMed]
24. EW, Ren W, Vanden-Eijnden E. Phys Rev B. 2002;66:52301.
25. Maragliano L, Fischer A, Vanden-Eijnden E, Ciccotti G. J Chem Phys. 2006;125:024106.
26. Pan AC, Sezer D, Roux B. J Phys Chem B. 2008;112:3432–3440. [PMC free article] [PubMed]
27. Huo S, Straub JE. J Chem Phys. 1997;107:5000–5006.
28. Vanden-Eijnden E, Venturoli M. J Chem Phys. 2009;130:194103. [PubMed]
29. Hummer G, Kevrekidis IG. J Chem Phys. 2003;118:10762–10773.
30. Berezhkovskii AM, Zitserman VY. Physica A. 1990;166:585.
31. Berezhkovskii AM, Berezhkovskii LM, Zitzerman VY. Chem Phys. 1989;130:55.
32. Best RB, Hummer G. Proc Natl Acad Sci USA. 2010;107:1088–1093. [PubMed]
33. Singer A, Erban R, Kevrekidis IG, Coifman RR. Proc Nat Acad Sci USA. 2009;106:16090–16095. [PubMed]
34. Branduardi D, Gervasio FL, Parrinello M. J Chem Phys. 2007;126:054103. [PubMed]
35. Ermak DL, McCammon JA. J Chem Phys. 1978;69:1352.
36. Vanden-Eijnden E. Lecture Notes in Physics. 2006;703:1–40.
37. EW, Vanden-Eijnden E. J Stat Phys. 2006;123:503–523.
38. Berezhkovskii A, Szabo A. J Chem Phys. 2005;122:014503. [PubMed]
39. Graham R. Z Physik B. 1977;26:397–405.
40. Risken H. The Fokker-Planck Equation. Springer; Berlin: 1996.
41. Peters B, Zimmermann NER, Beckham GT, Tester JW, Trout BL. J Am Chem Soc. 2008;130:17342–17350. [PubMed]
42. Schenter GK, Garrett BC, Truhlar DG. J Chem Phys. 2003;119:5828–5833.
43. Darve E, Pohorille A. J Chem Phys. 2001;115:9169.
44. Hänggi P, Talkner P, Borkovec M. Rev Mod Phys. 1990;62:251–341.
45. Tal FA, Vanden-Eijnden E. Nonlinearity. 2006;19:501–509.
46. Hummer G. New Jour Phys. 2005;7:34.
47. Im W, Roux B. J Mol Biol. 2002;319:1177–1197. [PubMed]
48. Knott BC, Duff N, Doherty MF, Peters B. J Chem Phys. 2009;131:224112. [PubMed]
49. Rhee YM, Pande VS. J Phys Chem B. 2005;109:6780–6786. [PubMed]
50. Zhu F, Hummer G. Proc Nat Acad Sci USA. 2010;107:19814–19819. [PubMed]
51. Moro GJ, Cardin F. Chem Phys. 1998;235:189–200.