Reaction Rates and Transition Pathways

213

rate. Minimizing this cost functional as a function of the path is approximately equivalent to minimizing the average time required to pass from reactant to product. A second important feature is that the temperature is included in the cost function. The reaction pathway found to be the pathway of maximum flux may depend on the temperature [35]. The temperature dependence may manifest itself simply as saddle point avoidance or more globally through a complete change in the dominant reaction mechanism.

This method has been applied to derive a multitude of paths for the coil-to-helix transition in polyalanine using a continuum solvation model [36].

C. Diffusional Paths

An alternative method for computing reaction pathways is based on the idea of diffusion [43] or a ‘‘noisy’’ dynamics [44,45]. A prescient review article by Wolynes [8] anticipated many of the most impressive advances in this area.

Consider a diffusion process in which a molecule is initially in a conformation defined by the coordinates rk. A transition probability can be constructed using a product of joint probabilities p(rk 1|rk) for moving between intermediate positions rk and rk 1 integrated over all intermediate states,

p(rP|rR) ∫drM 1 dr2 dr1 p(rP|rM 1) p(r2|r1)p(r1|rR)

(23)

Pratt [43] made the innovative suggestion that transition pathways could be determined by maximizing the cumulative transition probability connecting the known reactant and product states. That is, the most probable transition pathways would be expected to be those with the largest conditional probability.

In the presence of a potential U(r) the system will feel a force F(rk) rU(r)|rk. There will also be a stochastic or random force acting on the system. The magnitude of that stochastic force is related to the temperature, the mass of the system, and the diffusion constant D. For a short time, it is possible to write the probability that the system has moved to a new position rk 1 as being proportional to the Gaussian probability [43]

1	rk DδtβF(rk)]2
p(rk 1|rk) exp 4Dδt[rk 1		(24)

In the absence of an external force, the probability of moving to a new position is a spherically symmetrical Gaussian distribution (where we have assumed that the diffusion is spatially isotropic).

The diffusion constant should be small enough to damp out inertial motion. In the presence of a force the diffusion is biased in the direction of the force. When the friction constant is very high, the diffusion constant is very small and the force bias is attenuated— the motion of the system is strongly overdamped. The distance that a particle moves in a short time δt is proportional to √Dδt.

Recently, Chandler and coworkers [46,47] revisited this idea and developed an elegant and promising methodology for the computation of reaction pathways and transition rates in molecular systems.

D.Onsager–Machlup Paths

Olender and Elber [45] made a novel suggestion that dynamic trajectories for long time processes having known initial and final states may be computed by using ‘‘noisy’’ dy-

214

Straub

namic trajectories. They define a probability that a trajectory would move from position rk to rk 1 as

p(rk 1|rk) exp	1	[m(rk 1 2rk rk 1) δt2F(rk 1)]2	(25)
	2σ2

The highest probability paths will make the argument of the exponential small. That will be true for paths that follow Newtonian dynamics where mr¨ F(r). Olender and Elber [45] demonstrated how large values of the time step δt can be used in a way that projects out high frequency motions of the system and allows for the simulation of long-time molecular dynamics trajectories for macromolecular systems.

Elber et al. [48] applied this method to explore the dynamics of the C-peptide in water with impressive results. More than 30 trajectories of C-peptide were generated, and the process of helix formation in water was examined. Remarkably, a time step of 500 ps was used, which allowed for the study of peptide folding on extended time scales.

VI. HOW TO CONSTRUCT A REACTION PATH

In the computation of characteristic reaction pathways, it is essential to include a number of constraints and restraints on the reaction pathway. For example, for conformational transitions in a macromolecule, rigid-body translation and rotation should be constrained. In addition, to avoid clustering of intermediate configurations in potential energy minima, the distances between intermediate points along the path should be restrained to be roughly equal. Such a restraint forces the system to take steps of regular size along the reaction pathway.

In many cases, it is also helpful to have the path repel itself so that the transition pathway is self-avoiding. An actual dynamic trajectory may oscillate about a minimum energy configuration prior to an activated transition. In the computed restrained, selfavoiding path, there will be no clusters of intermediates isolated in potential energy minima and no loops or redundant segments. The self-avoidance restraint reduces the ‘‘wasted effort’’ in the search for a characteristic reaction pathway. The constraints and restraints are essential components of the computational protocol.

A. The Use of Constraints and Restraints

Following the computational protocol of Czerminski and Elber [39,40] a number of restraints and constraints are added to (1) encourage the mean-square distances between adjacent structures to be approximately constant,

M
CA(R) κ [(rk rk 1)2 dave2 ]2	(26)

k 1

where dave2 ∑Mk 1 (rk rk 1)2/M, (2) discourage intermediates from coming close to one another,

	ρ	exp	λ(rj rk)2
CR(R)	λ		d 2		(27)
		j k 1

Reaction Rates and Transition Pathways

215

and (3) eliminate rigid-body translations and rotations [39]

N
m (r rfix ) 0	(28)
1
and
N
m r rfix 0	(29)

1

where N is the number of atoms in the system, m is the atomic mass, and r gives the Cartesian coordinates for the µth atom. {rfix } 1,N is the arithmetic average of the coordinates of the atoms in the reactant and product configurations. In the final refinement, these terms do not add significantly to the integrated cost function. These constraints have been applied in a variety of contexts. Values for the various parameters can be found in a number of works [15,35].

B. Variationally Optimizing the Cost Function

The goal is to find the global minimum value of the reaction time function in the space of all possible reaction paths subject to the restraints. This is a computationally demanding task in dealing with a large biomolecule. One method is to minimize the function of reaction time using a conjugate gradient minimization of the chain of intermediates with the reactant (rR) and product (rP) configurations remaining fixed. Such a ‘‘focal method’’ leads to a solution that is strongly dependent on the goodness of the initial guess.

Alternative algorithms employ global optimization methods such as simulated annealing that can explore the set of all possible reaction pathways [35]. In the MaxFlux method it is helpful to vary the value of β (temperature) that appears in the differential cost function from an initially low β (high temperature), where the effective surface is smooth, to a high β (the reaction temperature of interest), where the reaction surface is more rugged.

Using a differential cost function such as that of Elber and Karplus, the potential energy is averaged over the path by including a factor of 1/L. In other definitions, such as the one employed in the MaxFlux method, there is no such normalization. Therefore, if the potential is set to zero, the MaxFlux method will find that the best path is the straight line path connecting reactants and products. However, methods where the differential cost is proportional to 1/L will find that all paths are equally good.

VII. FOCAL METHODS FOR REFINING TRANSITION STATES

Using the global search methods described above, it is possible to identify one or a number of important transition pathways. The coarseness of the search (the number of intermediate structures used to define the path) imposes a limit on smoothness of the pathway (see Fig. 8). In general, we can search along the reaction pathway and find the region of greatest energy curvature, | 2qE|max. There will also be an average distance between intermediates along the reaction coordinate δq2. The product of these numbers provides an estimate of the error in the energy of the pathway:

216

Straub

Figure 8 An accurate estimate of the barrier height can be found by adding a sufficient number of intermediate points in the discretized transition pathways. The solid line in the graph represents the energy profile for a reaction path described by 11 intermediate configurations of the system. The dashed line shows a coarse pathway described by only two intermediate configurations. The latter path underestimates the true energy barrier.

Maximum error | q2E| maxδq2

(30)

For a simple bistable reaction potential, it is clear that maximum curvature along the reaction pathway will occur near the extrema—the minima and the barrier top. The path endpoints are typically chosen to sit in the reactant and product minima, and in such a case the maximum error will result from the path ‘‘straddling’’ the barrier top as in Figure 8. Of course, this is the error made in a single segment of the pathway. For a general potential the pathway will consist of multiple segments and may have many barriers.

The curvature along the reaction coordinate is fixed by the energy function. It is possible to reduce the error by adding intermediate structures and limiting the error by reducing δq2, but there is an associated computational cost. It may be wiser to employ a focal method to refine the transition pathway—a local search for transition states and minima in the neighborhood of the globally defined pathway. A good number of very fine