Reaction Rates and Transition Pathways

209

tial energy is a symmetrical double well, P P P, and this result reduces to κ P/(2 P) as we found earlier.

IV. FINDING GOOD REACTION COORDINATES

In simple chemical systems, it is often possible to make a good first guess at the dominant reaction pathway [25–28]. An example of such a reaction is the chair-to-boat isomerization in cyclohexane. In that pathway, a clever combination of two torsion angles provides an excellent reaction coordinate for the isomerization reaction [29,30].

In more complex systems, the dominant reaction pathway may be less than obvious. This may be the case even for a modest configurational rearrangement. A fine example is the transport of a Na ion through the gramicidin ion channel. Even for such a welldefined process, with a highly constrained pathway for the ion to follow, a variety of reaction coordinates can be found with significantly different energetic features [31,32]. Moreover, even a minimal energy reaction pathway may include motion of atoms of the ion channel and water molecules in addition to the ion itself. Such a transition pathway may be difficult to express in terms of local atomic coordinates or simply defined collective modes.

A second example is the isomerization of a tyrosine residue about the ξ2 torsion angle in the protein bovine pancreatic trypsin inhibitor [33]. Largely buried in the protein core, the tyrosine 35 ‘‘ring flip’’ is an infrequent event occurring on the time scale of seconds owing to an activation energy barrier estimated on the order of 13 kcal/mol. The first guess for a good reaction coordinate, the ξ2 torsion angle, is inadequate. Exploring along this poor reaction coordinate, one grossly underestimates the energy barrier. What is missing? A section of backbone adjacent to the tyrosine ring is displaced in the transition. In fact, the barrier to rotation is largely determined by the steric contact between the rotating phenyl ring and the adjacent amide group, which must be included in the definition of the reaction coordinate [34].

If the motion of this group is ignored in the definition of the reaction coordinate, a search along the ξ2 angle alone can provide a misleading estimate of the reaction energetics. This is indicated in Figure 6. The potential is essentially two-dimensional. If the reaction coordinate is described too simply (using a local Cartesian or internal coordinate rather than a globally optimized reaction pathway), the potential energy may indicate an artificially low energy barrier. The discrete representation of the reaction path in terms of a set of intermediate structures renders the ‘‘motion’’ along the reaction pathway essentially discontinuous. (A poor choice of the reaction coordinate q, one that might not be a single-valued function of the computed reaction pathway l(r), can also lead to essential discontinuities in the reaction coordinate—points where dq/dl 0.) In the region of the transition state the ‘‘true’’ reaction coordinate is practically orthogonal to the chosen reaction coordinate. A significant change in energy along the true reaction coordinate occurs (a large jump upward), with almost no change in the chosen reaction coordinate. Proper use of both coordinates in the definition of the reaction coordinate will lead to a continuous pathway that represents well the energetics of dynamic barrier crossing.

In other cases, it may be impossible to describe the kinetics properly using a single reaction path. A variety of pathways may contribute to the reaction kinetics. One or more paths may be dominant at low temperature, whereas other paths may be dominant at high temperatures. This results in a temperature-dependent reaction mechanism. In such situa-

210

Straub

Figure 6 A poor choice of reaction coordinate can lead to a poor estimate of the activation energy and related rate constant. Because of the discrete nature of the reaction pathway, it is possible to step over the barrier. This leads to an underestimate of the activation energy.

tions, the goal is to be able to characterize the dominant reaction pathway or pathways at a given temperature [35–37].

V.GLOBAL SEARCHES FOR IMPORTANT TRANSITION PATHWAYS

In a typical dynamic trajectory, the initial position is well controlled but the endpoint of the trajectory is unknown. For chemical reaction dynamics, we are interested in trajectories that link known initial (reactant) and final (product) states so that both the initial conditions and the final conditions of the trajectory are fixed.

Reaction Rates and Transition Pathways

211

A. Variational Methods for Computing Reaction Paths

A most successful paradigm for isolating reaction pathways in complex systems begins with a definition of the reaction pathway as a continuous line l(r) that connects known reactant rR and product rP configurations. We then define an integrated cost functional

rP
cost[l(r)] rR (r)dl(r)	(19)

as a function of the path l(r) leading from the reactant configuration rR to the product configuration rP. The goal is to minimize this functional in the space of all possible paths. That is, we want to find the extremum of this function by varying the path l(r). The term(r)dl(r) assigns a particular differential cost or penalty for motion about the position r over an increment dl(r). We usually think of the cost function as a positive real number.

B. Choice of a Differential Cost Function

How can the continuous transition pathway l(r) be represented by a computer for a complex molecular system? l(r) can be approximately represented as a set of configurations of the system {rk}.

An example is shown in Figure 7 for the case of the coil-to-helix transition. The endpoints of the calculation are an unstructured coil rR and helix rP. Intermediate peptide structures correspond to transition intermediates defining the pathway l(r).

Elber and Karplus [38] presented an effective set of numerical methods for computing the reaction paths based on this approximation. First the path is discretized—it is expressed as a chain of intermediate configurations of the system {rk}. The line integrals of Eq. (19) are then written as

M 1
cost[l(r)] (rk)|rk 1 rk|	(20)

k 0

This discretized path allows us to represent the transition pathway numerically as a set of discrete configurations of the system.

1. Elber and Karplus Reaction Paths

What is the best choice of differential cost function? A variety of definitions of the cost function have been proposed. One stems form the highly original work of Elber and Karplus [38] and Czerminski and Elber [39], where

	(r)	1	U(r)	(21)
		L

U(r) is the potential energy of the molecular system, and L is the total length of the path L ∫dl(r). In this case, the cost is the average potential energy along the pathway. The best path is defined as the path of minimum mean potential energy. This definition has been applied to compute the reaction paths of a large number of molecular systems including peptide conformational rearrangement [40] and ligand migration in heme proteins. Similar methods have been shown to have greater numerical efficiency in applications to large systems [41].

212

Straub

Figure 7 The discretized reaction path is represented by a number of intermediate configurations of the system connecting the fixed reactant (1, coil) and product (12, helix) states.

2. MaxFlux Reaction Paths

An alternative definition of the cost function was proposed by Huo and Straub [35] based on an earlier suggestion of Berkowitz et al. [42],

C(r) eβU(r)

(22)

This result comes from the idea of a variational rate theory for a diffusive dynamics. If the dynamics of the reactive system is overdamped and the effective friction is spatially isotropic, the time required to pass from the reactant to the product state is expected to be proportional to the integral over the path of the inverse Boltzmann probability.

An important feature of this method is the connection to variational rate theory: The best reaction pathway is one that minimizes the reaction time or maximizes the reaction