Reaction coordinates and rates from transition paths

Abstract

The molecular mechanism of a reaction in solution is reflected in its transition-state ensemble and transition paths. We use a Bayesian formula relating the equilibrium and transition-path ensembles to identify transition states, rank reaction coordinates, and estimate rate coefficients. We also introduce a variational procedure to optimize reaction coordinates. The theory is illustrated with applications to protein folding and the dipole reorientation of an ordered water chain inside a carbon nanotube. To describe the folding of a simple model of a three-helix bundle protein, we variationally optimize the weights of a projection onto the matrix of native and nonnative amino acid contacts. The resulting one-dimensional reaction coordinate captures the folding transition state, with formation and packing of helix 2 and 3 constituting the bottleneck for folding.

Fig. 1.

Examples of “poor” and “good” reaction coordinates for the folding of a 47-residue Gō-like protein model. The poor coordinate (r_uni) is the projection of the contact map onto a matrix with uniform weights, whereas the good coordinate (r_opt) projects the contact map onto a matrix with optimized weights. (A and D) Time series of each reaction coordinate for the same simulation segment. (B and E) The equilibrium probability of the reaction coordinate p_eq(r). (C and F) The probability of being on a transition path given the value of r, p(TP|r). The green horizontal lines indicate “transition states” r = r^‡. Only 1/10th of the trajectories were used in the optimization to avoid “over-fitting,” but the complete trajectories were used to test r with p(TP|r). Insets in C and F show the distributions of splitting probabilities for configurations at the maximum r = r^‡ of p(TP|r).

Fig. 2.

Structures representing different values of the folding coordinate r_opt. Residues are colored on a blue–green–red color scale from the N to C terminus and aligned with the native structure by using residues 17–47 (helices 2 and 3) only. (A) Folded side of the barrier (r_opt = 7.4). (B) Transition state at the maximum of p(TP|r_opt) (r_opt = 6.9). Note that helix 1 is at least partially formed, but not properly docked against the scaffold of helices 2 and 3. (C) Unfolded side of the barrier (r_opt = 6.5). (D–F) Fraction of contacts 〈q_ij〉 present in the folded state (r_opt > 7.2) (D), transition state (6.89 < r_opt < 6.91) (E), and unfolded state (r_opt < 6.5) (F). Note that helical contacts are overemphasized by the 12-Å contact cutoff.

Fig. 3.

Dipole flip of water chain inside carbon nanotube. (A) Equilibrium free energy surface as a function of the total dipole moment (D, debye unit) of the water chain from umbrella sampling (red line) and equilibrium MD (blue squares). (B) Probability density of the dipole moment in the transition-path ensemble. (C) Transition-path probability p(TP|M_z). (D) Reciprocal rate coefficient k^–1 for dipole flip from ≈66-ns equilibrium MD simulations (gray shaded area indicates ± one estimated standard deviation), 〈t_TP〉/p(TP) (blue rectangle indicates ± one estimated standard deviation; red symbols with error bars show estimates from individual histogram bins), and transition-state theory (green line) (3, 42).

Fig. 4.

Snapshots of water-chain reorientation inside carbon nanotube along the transition path with the highest relative weight, Eq. 5. Blue arrows indicate the progression of a hydrogen-bond defect along the water chain over a period of 0.7 ps.