2.1 A simple model with intra-state equilibration and inter-state jump

As noted above, a state consists of a set of local energy minima. Each local minimum corresponds to a microstate. Figure 2a presents a model in which state A consists of three microstates, to be referred to by index i or j, with i and j=1–3; from microstate 3 the molecule can make a jump to state B. We choose model parameters to illustrate two situations: and . The results of these two situations are displayed in Fig. 2b and c, respectively.

Fig. 2.

A simple model with intra-state equilibration and inter-state jump. (a) Illustration of the model. The two states, A and B, are represented by boxes; microstates within state A are represented by circles. (b) Time dependences of σ₁ and ρ₁/ρ_Aeq, representing intra-state equilibration, and of ρ_A, representing inter-state jump, for the following parameter values: k ₁₂=k ₂₁=k ₁₃=k ₂₃=1, k ₃₁=k ₃₂=2, and k _tr=0·1. (c) The corresponding results when k ₃₁ and k ₃₂ are decreased to 0·125 and k _tr is increased to 1. The results are obtained from kinetic simulations of a single molecule. Briefly, the waiting time of the molecule in an initial microstate is generated from an exponential distribution function [see Eq. (2.10)], with the average waiting time equal to the inverse of the sum of the rate constants for all the pathways leaving that microstate. The probability for taking each of these pathways is proportional to the corresponding rate constant. The results shown are the average of 10⁶ repeat simulations; each simulation starts with the molecule in microstate 1. This simulation procedure is similar in spirit to the stochastic simulation algorithm of Gillespie (Reference Gillespie1977).

First, let us study the equilibration among the three microstates of state A, which we model as rate processes, with rate constants k_ij. For now jumps to state B are not considered. The rate constants for the forward and reverse transitions between two microstates, e.g., k ₁₂ and k ₂₁, are constrained by the corresponding equilibrium constant, which in turn is given by the ratio of the configurational integrals of the two microstates [Eq. (2.11)]. Our interest is in how the proportions of molecules in the microstates, hereafter referred to as occupation probabilities, evolve over time. We denote the occupation probabilities here as σ_i; the corresponding quantities when jumps to state B are allowed will be denoted as ρ_i (see below). Starting from any initial values σ_i0, σ_i will relax to their equilibrium values σ_ieq, which are proportional to the respective configurational integrals. The time dependence of σ_i can be obtained from a kinetic simulation (see Fig. 2). In Fig. 2b, we plot σ₁−σ_1eq as a function of time when the initial values are σ₁₀=1 and σ₂₀=σ₃₀=0 and the rate constants are k ₁₂=k ₂₁=k ₁₃=k ₂₃=1 and k ₃₁=k ₃₂=2. The corresponding equilibrium values are σ_1eq=σ_2eq=0·4 and σ_3eq=0·2. The equilibration time τ_Aeq can be estimated as the area under the (σ₁−σ_1eq)/(σ₁₀−σ_1eq) versus t curve; we find τ_Aeq=0·31. Figure 2c displays the result for σ₁−σ_1eq when k ₃₁ and k ₃₂ are decreased by 16-fold to 0·125 (while all other parameters are unchanged). The equilibrium values become σ_1eq=σ_2eq=0·2 and σ_3eq=0·8, and τ_Aeq increases to 0·54.

We now include jumps to state B from microstate 3, again modeled as a rate process, with rate constant k _tr. Both the equilibration among the three microstates and the jumps to state B affect ρ_i, the occupation probabilities. The latter leads to a decrease in the total occupation probability of state A, ρ_A=ρ₁+ρ₂+ρ₃, over time. The decrease is governed by

(2.12)

We use the normalized occupation probabilities, ρ_i/ρ_A, to monitor the equilibration among the three microstates. The situation displayed in Fig. 2b has k _tr=0·1, along with the values of k_ij specified above. Firstly, it can be seen that the decay of ρ_A can be fitted well to an exponential function, which is expected if a rate description is valid [Eq. (2.4)]. The fitted rate constant, k ₊, has a value of 0·0196; correspondingly, the average waiting time is 51. Secondly, the time dependence of ρ₁/ρ_A matches that of σ₁. Thus, the equilibration time of 0·31 based on σ₁ is also correct for the present case where jumps to state B are allowed. Compared to , we find that τ_Aeq is two orders of magnitude shorter in the present situation. That the initial decays of ρ_i/ρ_A and σ_i match are easy to rationalize, since during this period very few jumps to state B occur. On the other hand, that ρ_i/ρ_A settle to the same equilibrium values σ_ieq suggests a scenario that, when the equilibrium within state A is perturbed by a jump from microstate 3 to state B, it is recovered by re-equilibration among the microstates before the next jump occurs. This scenario also allows us to find an approximate value for the rate constant k ₊. Making use of the result ρ₃/ρ_A≈σ_3eq, which holds after a brief initial period, Eq. (2.12) becomes

(2.13a)

The rate constant is thus

(2.13b)

which has a value of 0·02 for the parameters used in Fig. 2b. This predicted value is very close to the value, 0·0196, obtained above from fitting the time dependence of ρ_A to an exponential function.

The situation displayed in Fig. 2c has k _tr=1. Compared to the parameters used in Fig. 2b, k _tr is increased by 10-fold, along with the 16-fold decrease in k ₃₁ and k ₃₂ stated above. For these parameters, the time dependence of ρ_A no longer fits well to an exponential function, and hence a rate description is not valid here. At the same time, ρ₁/ρ_A deviates significantly from σ₁, except at short times (when jumps to state B are only few). Notably, ρ₁/ρ_A and ρ_A decay on a comparable timescale, validating the claim that the case under consideration has .

Together, the contrasting situations shown in Fig. 2b, c illustrate the point that a rate description will hold when the transitions between two states are rare enough to allow for intra-state equilibration.

2.2 Intermolecular transitions

Up to now, we have discussed transitions between states of the same molecule. Many reactions of biological interest involve the binding of two different molecules, e.g., a protein and a ligand, to form a complex. In this case, the protein molecule makes a transition from a state in which it is freely translating and rotating by itself throughout the container to a state in which its translation and rotation are coupled to those of the ligand molecule. (The alternative perspective, centered on the ligand instead of the protein, is equally valid.) Even though the transition is intermolecular instead of intramolecular, the derivation leading to Eq. (2.2) can still be considered valid, except that the ‘rate constant’ k introduced there is now proportional to the ligand concentration C _L, since the probability that a protein molecule makes the transition from the free state to the protein–ligand complex in an infinitesimal time interval is proportional to C _L. We thus replace k by k _aC _L, where k _a is the new rate constant for the association of the protein molecule with the ligand molecule. The reverse process, i.e., the dissociation of the protein–ligand complex into separate protein and ligand molecules, is an example of the intramolecular transitions discussed above, since the reactant, i.e., the protein–ligand complex, is effectively one molecule. We denote the dissociation rate constant as k _d. The association and dissociation can be represented by the reaction scheme

(2.14)

The rate equation for the protein concentration, C _P, is

(2.15)

Here the forward reaction rate is second order in reactant concentration and hence that reaction is referred to as second order. When C _L is a constant, which, e.g., is practically realized when the ligand is in large excess over the protein, the product k _aC _L can be treated as an effective rate constant. For such a pseudo-first-order situation, the solution of the rate equation can be obtained by identifying k _± in Eq. (2.8a) with k _aC _L and k _d, respectively, leading to

(2.16)

where we have introduced the equilibrium constant for dissociation, or dissociation constant

(2.17)

Its inverse is called the association constant, denoted as K _a. The equilibrium protein concentration is given by

(2.18)

where C _Pt is the total protein concentration, either in the free state or in the bound state:

(2.19a)

and ρ_beq gives the bound fraction (the unbound fraction ρ_ueq is 1−ρ_beq).

When C _L cannot be treated as a constant, it is still constrained by the fact that the total ligand concentration, either free or in the complex, is fixed:

(2.19b)

Using Eqs. (2.19) in Eq. (2.15), we find the solution for C _P to be given by

(2.20)

where

(2.21a)

(2.21b)

It can be verified that when C _Pt≪C _Lt, Eq. (2.20) reduces to Eq. (2.16).

At long times, the equilibrium value that the protein concentration settles into is C _P−. Unlike the unimolecular case [Eq. (2.8)], where the equilibrium concentration as a proportion of the total concentration is an intrinsic property of the molecule, in the bimolecular case the equilibrium protein concentration as a proportion of the total protein concentration is not an intrinsic constant; it depends on the total protein concentration and the total ligand concentration. Similar to Eq. (2.11) for the unimolecular case, we may express the association constant in terms of the configurational integrals of the protein, the ligand, and their complex. Let the potential energy functions of the three systems be U _P(x_P), U _L(x_L), and U _PL(x_PL). Then we have (Zhou & Gilson, Reference Zhou and Gilson2009)

(2.22)

In some reactions, product formation involves more than two reactant molecules. For example, many enzymatic reactions require the binding of two or more substrates to one enzyme. Compared to the binding of two molecules to form a binary complex, the chance of three molecules binding simultaneously to form a ternary complex is negligibly small. Surely the ternary complex must form in two steps: in the first step, two of the molecules bind to form an intermediate binary complex; in the second step, the third molecule binds to the pre-formed binary complex to yield the ternary complex. It can thus be concluded that all elementary reactions are either unimolecular or bimolecular.

3.1 Eyring's transition-state theory

The most influential rate theory is the transition-state theory of Eyring (Reference Eyring1935), also worked out by Evans & Polanyi (Reference Evans and Polanyi1935). It is based on three assumptions:

1. (i) The motion of molecules obeys Newton's equation. The potential energy function for our illustrative model is shown in Fig. 3a.

2. (ii) As soon as a molecule crosses the energy barrier at x=x‡, or transition state, it is considered to form the product.

3. (iii) Molecules in the barrier region are in thermal equilibrium with those in the energy well of the reactant state.

In the barrier region, the potential energy function is flat; hence, there the force is zero and molecules undergo free translation. If δN molecules, all with velocity v>0, are found in the region with length δl just to the left of x=x‡, then within the time interval dt each molecule moves a distance of vdt to the right and the number of molecules that cross the barrier is δN⋅(vdt/δl). Consequently, the rate constant is

(3.6)

where N is the number of molecules found in the reactant energy well.

Under thermal equilibrium, δN/N is the same as the ratio of the partition functions in the barrier region and in the reactant energy well:

(3.7)

The partition function in the barrier region is that appropriate for free translation in one dimension:

(3.8)

where p=mv is the momentum of a molecule along the x-direction, m is the mass, and ΔU‡ is the activation energy (Fig. 3a). Following convention we have inserted Planck's constant h to make the partition function unitless. As indicated by Eq. (3.8), the molecules in the barrier region have a continuous distribution of velocities; therefore, the velocity v appearing in Eq. (3.6) should be replaced by its average:

(3.9)

where θ(v), a Heaviside function with value 1 if v>0 and 0 otherwise, is used to select only positive velocities. Combining the above results, we find

(3.10)

We now discuss the calculation of q _A. If the problem is treated according to classical statistical mechanics, we have

(3.11)

Equation (3.10) becomes

(3.12a)

Note that Planck's constant drops out of the final result for k ₊, as to be expected from a classical treatment. For later use, we rewrite it in the form

(3.12b)

where 〈···〉 denotes an average with the canonical-ensemble equilibrium probability density

(3.13)

δ(x) is a delta function, and

(3.14)

is the equilibrium occupation probability in state A. The corresponding result for k ₋ is obtained by replacing ρ_Aeq with ρ_Beq in Eq. (3.12b). It can then be easily seen that the ratio k ₊/k ₋ is the same as the equilibrium constant ρ_Beq/ρ_Aeq, as demanded by Eq. (2.9). We may further make a harmonic approximation of the potential energy function around the minimum x _A of the reactant well:

(3.15)

Using this approximation and evaluating the configurational integral in Eq. (3.12a), the rate constant is now given by

(3.16)

where ω_A=f _A^1/2 is the angular frequency of the harmonic oscillation around x=x _A.

With the harmonic approximation of the potential energy function, a quantum treatment for calculating q _A is possible. Essentially q _A is the partition function of a harmonic oscillator, given by

(3.17)

where ħ=h/2π. Equation (3.10) is now

(3.18)

This result was first derived by Herzfeld (Reference Herzfeld1919) for the dissociation of a diatomic molecule. In the high-temperature limit, q _A→k _BT/ħω_A, and we recover the classical result of Eq. (3.16). This is to be expected, since quantum statistical mechanics approaches the classical limit at high temperatures.

Below we further comment on the approximation of the transition-state theory and a common misunderstanding of the theory.

3.1.1 Approximation of the transition-state theory

The nature of the approximation of the transition-state theory can be elucidated in the context of a more exact microscopic formulation of the rate constant. Equation (2.8a), the solution of the rate equation, states that the relaxation of ρ_A, the occupation probability in the reactant state, toward its equilibrium value is an exponential, with the relaxation time τ_rxn=(k ₊+k ₋)⁻¹. Microscopically ρ_A is the average of a variable that has value 1 if the molecule occupies the reactant well and 0 otherwise:

(3.19)

Now Onsager's regression hypothesis states that the regression of microscopic thermal fluctuations at equilibrium follows the macroscopic law of relaxation of small non-equilibrium disturbances (Onsager, Reference Onsager1931). Accordingly,

(3.20)

where

(3.21)

is the correlation function of the instantaneous fluctuation θ(x‡−x(t))−ρ_Aeq, and is referred to as the number correlation function. Noting that θ²(x)=θ(x) and 1−ρ_Aeq=ρ_Beq, it can be easily shown that c(0)=ρ_Aeqρ_Beq. Now Eq. (2.8a) becomes

(3.22)

To make connection with the transition-state theory, we take the time derivative of Eq. (3.22):

(3.23)

where a dot represents a derivative. The time derivative of the number correlation function is given by

(3.24)

which can be interpreted as the reactive flux, i.e., the flux across the dividing surface given that the reactant molecule starts there and is found in the product well time t later (Chandler, Reference Chandler1978). As discussed in section 2, a rate description is valid only when there is a significant separation in timescales between intra-state equilibration (τ_eq) and inter-state transitions (τ_rxn). In that case, the number correlation function is exponential at times longer than τ_eq. Because of the separation in timescales, there exist intermediate times, denoted as Δt, such that τ_eq≪Δt≪τ_rxn. At these intermediate times, exp(−Δt/τ_rxn) is very close to 1; therefore Eq. (3.23) becomes (Chandler, Reference Chandler1978; Yamamoto, Reference Yamamoto1960)

(3.25a)

Correspondingly,

(3.25b)

The transition-state theory amounts to setting Δt in the last equation to 0⁺, which denotes a time that is greater than 0 by an infinitesimal amount:

(3.26)

The equivalence of Eqs. (3.26) and (3.12b) is seen when one realizes that x(0⁺)>x‡ if and only if v(0)>0. A related approach for formulating the rate constant is presented in subsection 3.2.

To recapitulate, the transition-state theory predicts the rate constant from the initial reactive flux, when the exponential relaxation of the number correlation function has not been established. The correct formulation is to use the reactive flux at the intermediate times Δt. Since the number correlation function decays with time, it can be expected that the transition-state theory always overestimates the rate constant. Although Eyring recognized the overestimation and introduced a transmission coefficient, it was not possible to calculate the transmission coefficient within the transition-state theory itself. That the transition-state theory provides an upper bound for the rate constant is the basis of variational transition-state theory, which seeks a dividing surface that minimizes the rate constant calculated by the transition-state theory (Keck, Reference Keck1960; Wigner, Reference Wigner1937).

We end this subsubsection with some results that are related to the preceding formulation and will become useful below. When the motion of the reaction coordinate is governed by the Langevin equation, Onsager's regression hypothesis [Eq. (3.20)] has been proven explicitly (Zhou, Reference Zhou1989). It is assumed the initial probability density is an equilibrium distribution confined to the reactant state. Then

(3.27a)

(3.27b)

where ρ_N(x,v,t|x ₀,v ₀) is the conditional probability density for finding the molecule at x with velocity v at time t given that it was at x ₀ with velocity v ₀ at t=0; ρ_N(x,t|x ₀) is the reduced conditional probability density in x for an equilibrium distribution in v ₀. The subscript ‘N’ has a special meaning: it signifies that the probability density satisfies natural boundary conditions, i.e., the decay to zero at x=±∞ because of the infinite potential energies there. These boundary conditions are to be distinguished from, e.g., absorbing boundary conditions that are introduced below for various purposes. By design ρ_A(0)=1. Instead of finding k ₊+k ₋ from the time derivative of the number correlation function, its inverse, i.e., the relaxation time, is given by the time integral [see Eq. (3.22)] (Skinner & Wolynes, Reference Skinner and Wolynes1978):

(3.28)

Now expressing c(t) in terms of ρ_A(t) [Eq. (3.20)], we find

(3.29)

Correspondingly (Zhou, Reference Zhou2005a),

(3.30)

3.1.2 Multi-dimensional transition-state theory and a common misunderstanding

Although we have so far presented the transition-state theory for a one-dimensional potential energy surface, the theory was originally derived for molecules described by a multi-dimensional energy function, to be denoted as U(x). It is assumed that a suitable reaction coordinate, x, is found; a surface, x=x‡, is appropriate as a dividing surface between the reactant and the product states, and x, at least around the minimum of U(x) on the dividing surface, is uncoupled to the remaining degrees of freedom. Now the whole dividing surface x=x‡ is the transition state. Let the partition function in the reactant well, with its minimum chosen as the reference of potential energy, still be denoted as q _A and the partition function for all coordinates other than x, with x fixed at x‡ and the minimum of U(x) on the dividing surface chosen as the reference of potential energy, be denoted as q′‡. Then Eq. (3.10) for the rate constant is generalized to

(3.31)

where the activation energy ΔU‡ is the difference between the minima of U(x) on the dividing surface and in the reactant well. The prime in the partition function q′‡ is very important; again, it signifies that the reaction coordinate is left out of the calculation of this partition. In the one-dimensional model discussed above, there are no coordinates other than the reaction coordinate, and q′‡ would be 1.

Often q′‡ is mistakenly thought as the partition function of the transition state. Following this mistaken notion, one would be misled to think that defines the free-energy difference, to be denoted as ΔG‡, between the transition state and the reactant state:

which in turn transforms Eq. (3.31) into

While mathematically the last expression cannot be faulted, it nevertheless results in the misleading interpretation that ΔG‡ is the activation energy and k _BT/h is the pre-exponential factor. The presence of Planck's constant in such a pre-exponential factor would suggest that quantum effects persist at all temperatures. In reality, q _A involves one additional coordinate than q′‡; at high temperatures this additional coordinate would contribute a factor that cancels the h in k _BT/h, as illustrated above on the one-dimensional model.

We again stress that in calculating q′‡ the reaction coordinate is left out, and q′‡ is not the partition function of the transition state. In fact, for the one-dimensional model, we have defined a partition function of the transition state [Eq. (3.8)], with the reaction coordinate not fixed at a single value (as would be required in calculating q′‡) but extending over a range δl. A rigorous implementation of the multi-dimensional transition-state theory will be illustrated below in subsection 3.3.

3.2 Kramers' turnover problem

In a sense, assumptions (i) and (iii) that lead to the transition-state theory are contradictory to each other. If the motion of molecules, each with a single degree of freedom along the reaction coordinate, is governed by Newton's equation, then the energy of each molecule would be conserved. Low-energy molecules would stay in the reactant well, and higher-energy molecules would stay in the barrier region; these two types of molecules would not be able to mix and equilibrate. In reality, the reaction coordinate is coupled to other degrees of freedom of a reactant molecule and to solvent molecules, which effectively serve as a thermal bath that allows for energy dissipation and thermal equilibration. In Kramers' theory (Kramers, Reference Kramers1940), the motion along the reaction coordinate is governed by the Langevin equation:

(3.32)

which augments Newton's equation by two additional forces: a frictional force, −mγv, which serves to dissipate the energy of the molecule,; and a random force, , which serves to agitate the molecule, leading to thermal equilibration. The random force has a Gaussian distribution, with zero mean and correlation given by

(3.33)

which is a form of the fluctuation–dissipation theorem.

Because of the random force, the trajectory of the molecule is not deterministic but stochastic. Complementary to stochastic trajectory, we may describe the microscopic behavior of the molecule by the probability density, ρ(x,v,t), in position and velocity as a function of time. This probability density satisfies the Fokker–Planck equation

(3.34)

It can be checked that the canonical-ensemble equilibrium distribution given by Eq. (3.13) is a stationary solution of the Fokker–Planck equation.

Let us first discuss the limits of high and low frictions. When γ→∞, the change in x is slow and the relaxation in v is fast. Then the inertial term on the left-hand side of the Langevin equation can be neglected, leading to

(3.35)

which is the governing equation for Brownian (or diffusive) motion of x.

Correspondingly, the probability density in x, ρ(x,t), satisfies the Smoluchowski equation

(3.36)

where D=k _BT/mγ is the positional diffusion constant. The equilibrium distribution is

(3.37)

In the high-friction limit, positional diffusion from the reactant well to the product well becomes rate limiting. It can thus be expected that in this limit k ₊∝D∝1/γ.

In the low-friction limit, the total energy, E=mv ²/2+U(x), is nearly conserved, and the motion of the molecule can be described as diffusion among different energy levels. The probability density in E, ρ(E,t), satisfies the equation (Zwanzig, Reference Zwanzig2001)

(3.38)

where

(3.39a)

is the action for the periodic orbit with total energy E and

(3.39b)

is the angular frequency of the orbit. Note that I(E) is also the area inside the orbit in phase (i.e., x and p) space, and dI(E)/dE=2π/ω(E) is the density of states. Here the equilibrium distribution ρ_eq(E)∝[2π/ω(E)]exp(−E/k _BT). Comparing Eqs. (3.38) and (3.36), it can be seen that the energy diffusion constant, D_E, is proportional to γ. Therefore, in the low-friction limit one expects k ₊∝D_E∝γ.

Putting the results in the two limits together, one sees that the rate constant increases with increasing γ when γ is small and decreases with increasing γ when γ is large (Fig. 3b). It can thus be expected that the value of k ₊ reaches a maximum at an intermediate γ and decreases to zero when γ approaches either zero or infinity. This dependence of the rate constant on friction is known as Kramers' turnover.

Kramers derived expressions for the rate constant for different ranges of the friction coefficient. His approach can be explained from the perspective of a hypothetical bulk experiment (Farkas, Reference Farkas1927; Zhou, Reference Zhou2005a). In this experiment, molecules are prepared at the reactant state with a concentration C. Whenever a transition is made to the product state, the molecule is immediately removed, as if there is an absorbing trap, and the reactant pool is replenished so that the reactant concentration is always maintained at C. After a brief transient period (lasting ~τ_eq), a steady state will set in. The steady-state flux will be [see Eq. (2.2)]

(3.40a)

Rearranging, we find

(3.40b)

The rate constant is thus given by the ratio of the steady-state flux and the reactant population. Note that the reactive-flux formulation, Eq. (3.25b), of the rate constant can be viewed as a special case of the flux-over-population formulation. As observed already, the transition-state theory [Eq. (3.26)] makes the approximation of replacing the ‘steady-state’ reactive flux by its initial value.

To implement the flux-over-population formulation, one would find a steady-state probability density that maintains the equilibrium distribution in the reactant well and satisfies the absorbing boundary condition in the product well. The implementation is straightforward in the high-friction regime, where the motion of the molecule is modeled as positional diffusion. A steady-state probability density corresponds to a constant flux:

(3.41)

The appropriate boundary conditions are

(3.42a)

(3.42b)

Solving Eq. (3.41) for ρ_ss(x) subject to the condition of Eq. (3.42a), we find

(3.43a)

Using Eq. (3.42b), we find the flux

(3.43b)

Now with the condition in Eq. (3.42a), the population in the reactant well is

The rate constant is finally

(3.44)

Making the harmonic approximation of Eq. (3.15) for the potential around x _A and a similar approximation

(3.45)

for the potential around the top of the barrier, we find

(3.46)

where ω‡=(f‡)^1/2. Equation (3.44) is still valid when D is a position-dependent diffusion coefficient, i.e., D(x). The first integral in the denominator is then dominated by the region around the maximum of an effective potential, U(x)−k _BTlnD(x). The location of the maximum, which was referred to as the kinetic transition state (Chahine et al. Reference Chahine, Oliveira, Leite and Wang2007), and the effective activation energy can be quite different from those of the original potential.

As shown below in subsubsection 3.6.2, the ratio ρ_ss(x)/ρ_eq(x) is the splitting probability, denoted as ζ_A(x), which is the probability that a molecule, started at x, will first reach the reactant well bottom instead of the product well bottom. When x is in the reactant well, ζ_A(x) is very close to 1. Around x=x‡, ζ_A(x) undergoes a rapid transition, to ~1/2 at x=x‡, and to 0 a small distance thereafter.

By a similar procedure, we find the rate constant in the low-friction limit:

(3.47a)

(3.47b)

where

(3.48)

with x _tA being the left-hand-side turning point, i.e., the solution of U(x)=ΔU‡ to the left of x=x _A. A comment on the upper limit of the second integral in Eq. (3.47a) is in order. Orbits with energies higher than ΔU‡ will traverse both the reactant and the product wells. Setting the upper limit of the integration to ΔU‡ amounts to assuming that all molecules moving along those orbits would be trapped in the product well. However, it should be recalled that, in the flux-over-population formulation, the absorbing boundary is set around x=x _B, i.e., inside the product well. For molecules with E>ΔU‡, this absorbing boundary traps only those moving to the left (i.e., v<0) but not those moving to the right (i.e., v>0). This ‘split’ boundary condition is treated properly later by Melnikov & Meshkov (Reference Melnikov and Meshkov1986); see subsection 3.4.

For moderate to high frictions, we again seek a steady-state probability density ρ_ss(x,v), which takes the equilibrium distribution around the bottom of the reactant well and becomes zero at the bottom of the product well.Footnote ¹ The flux from the reactant well to the product well is

(3.49)

We write ρ_ss(x,v) as

(3.50)

The splitting probability ζ_A(x,v)=ζ(x,–v) (the appearance of the negative sign will be explained in subsubsection 3.6.2). ζ(x,v) satisfies the equation

(3.51a)

with the boundary conditions

We now solve for ζ(x,v) around x=x‡. Applying the harmonic approximation of Eq. (3.45), Eq. (3.51a) becomes

(3.51b)

The trick is to seek a set of parallel lines in phase space, over each of which ζ(x,v) is a constant; the value of this constant changes as the line moves from the reactant well to the product well. Let each line be represented by the equation λ(x−x‡)−v=u, where u dictates the location of the line. The desired condition is satisfied by choosing

(3.52)

Correspondingly, the solution of Eq. (3.51b) is

(3.53)

which has a transition in value around u=0, i.e., the line that goes through the point (x, v)=(x‡, 0). If this transition is essentially complete within the range where the harmonic approximation of the potential around x=x‡ holds, then the boundary conditions are satisfied by setting

The rate constant can finally be calculated as

(3.54)

It can be easily checked that Eq. (3.54) reduces to the high-friction limit given by Eq. (3.46) when γ→∞. However, when γ→0, Eq. (3.54) does not reduce to the low-friction limit given by Eq. (3.47b). Instead, in this limit it reduces to the transition-state theory result given by Eq. (3.16). The validity of Eq. (3.54) requires that the transition of ζ_A(x,v) occurs in a range where the harmonic approximation of the potential around x=x‡ holds. As noted by Melnikov & Meshkov (Reference Melnikov and Meshkov1986), this range is proportional to γ^−1/2, which extends to infinity as γ→0. That explains why Eq. (3.54) is only valid for moderate to high frictions (see also footnote ¹).

In summary, Kramers has shown that the rate constant is proportional to γ when γ is low and proportional to γ⁻¹ when γ is high. The rate constant exhibits a maximum at an intermediate friction. However, Kramers did not find an expression for the rate constant that works for the full range of friction.

3.3 Friction with memory: Grote–Hynes theory

The Langevin equation (3.32) is a special case of a more general class of equation of motion

(3.55)

known as the generalized Langevin equation. The friction coefficient here has the form of a memory kernel. Correspondingly, the correlation function of the Gaussian random force is

(3.56)

Grote & Hynes (Reference Grote and Hynes1980) used the generalized Langevin equation to model the motion along the reaction coordinate and derived the counterpart to Kramers' result [Eq. (3.54)] for the rate constant at moderate to high frictions.

The Grote–Hynes theory is based on a formulation of the rate constant derived from the so-called stable states picture (Northrup & Hynes, Reference Northrup and Hynes1980). As explained in the preceding subsection, the steady-state probability density appropriate for calculating the rate constant is nearly identical to the equilibrium distribution in the reactant well but deviates significantly from it in the barrier region. Northrup & Hynes (Reference Northrup and Hynes1980) thus introduced an intermediate state, I, in addition to the ‘stable’ states A and B. The dividing line between A and I is located at a position where the steady-state probability density starts to deviate from the equilibrium distribution; a similar dividing line can be defined between B and I. These two positions will be denoted as x _A−I and x _I−B, respectively. The rate constant is given by

(3.57a)

where the first flux is across the A–I dividing line at time 0 and the second flux is across the I–B dividing line at time t. The initial probability density is an equilibrium distribution in A; an absorbing boundary is set in B. Equation (3.57a) is closely related to the reactive-flux formulation [Eq. (3.25b)] of the rate constant.

Similar to Kramers, Grote and Hynes assumed that the harmonic approximation of Eq. (3.45) applies to the entire region of I. They then extended this downward parabola indefinitely beyond x=x _A-I and x=x _I-B. Let the conditional probability density in this extended parabola be ρ_N(x,v,t|x ₀,v ₀). Note that ρ_N(x,v,t|x ₀,v ₀) naturally satisfies the desired absorbing boundary condition. The rate constant can now be expressed as

(3.57b)

The conditional probability density in a parabolic potential is Gaussian. Evaluating all the integrals involved, the final result is

(3.57c)

where the Grote–Hynes ‘reactive’ frequency ω_GH is the positive root of the equation

(3.58)

with denoting the Laplace transform of the memory kernel. For the memory-free friction, ξ(t)=2γδ(t), Eq. (3.58) gives , and one recovers Kramers' result, Eq. (3.54).

As alluded to in the preceding subsection, friction arises from coupling of the reaction coordinate to other degrees of freedom of the reactant molecule and to solvent molecules, serving as a thermal bath. A simple model for a thermal bath consists of a set of n harmonic oscillators, with the potential energy (Zwanzig, Reference Zwanzig2001)

(3.59)

The whole system, reaction coordinate plus thermal bath, satisfies Newton's equation. By eliminating the bath coordinates, y_j, one arrives at the generalized Langevin equation (3.55). The memory kernel is given by

(3.60)

and the ‘random’ force is given by

(3.61)

where x ₀ is the initial value of x, y_{i 0} are the initial values of the bath coordinates, and v_{i 0} are the corresponding initial velocities. Equation (3.60) shows that memory kernels, such as that of an exponential form, typically used to model unimolecular reactions correspond to an infinite set of bath oscillators that spans a continuous spectrum of frequencies. For later reference, we note that the Laplace transform of the memory kernel in Eq. (3.60) is given by

(3.62)

With a small set (n~30) of oscillators, computer simulations (Zhou & Zwanzig, Reference Zhou and Zwanzig2002) show that the Grote–Hynes theory provides a reasonable estimate for the rate constant at moderate frictions. However, at high frictions, the number correlation function becomes non-exponential, and the decay is much slower than that expected from the Grote–Hynes theory.

Interestingly, Pollak (Reference Pollak1986) found that, when the multi-dimensional transition-state theory [Eq. (3.31)] is applied to the whole system of reaction coordinate plus thermal bath, the Grote–Hynes result, Eq. (3.57c), is derived. We outline this derivation here, mainly to illustrate the use of the multi-dimensional transition-state theory. The implementation of Eq. (3.31) entails evaluating two partition functions, q _A and q′‡. To evaluate q _A, we again use the harmonic approximation of Eq. (3.15). Then the system is equivalent to a set of n+1 harmonic oscillators. The eigenvalues, λ_Ai, i=0 to n, of the [mass-scaled; cf. Eq. (3.81)] force-constant matrix, f_A, give the angular frequencies as λ_Ai^1/2. In the classical limit, each oscillator contributes a factor k _BT/ħ_Ai^1/2 to the partition function. Therefore,

(3.63)

where det denotes determinant. To calculate q′‡, we apply the harmonic approximation of Eq. (3.45) and carry out a similar normal mode analysis around the saddle point of the multi-dimensional energy function U(x)+U _bath(y). The force-constant matrix f‡ in the barrier region has one negative eigenvalue, −λ₀‡, and n positive eigenvalues, λ_i‡, i=1 to n. The value of λ₀‡ is determined by

(3.64)

where I denotes the identity matrix. The last equation is satisfied when

(3.65)

Comparison with Eq. (3.58) shows that =ω_GH. We define the mode corresponding to the negative eigenvalue as the modified reaction coordinate. This mode, with coordinate denoted as x′, does not coincide with the original reaction coordinate x and instead is a linear combination of x and the bath coordinates. For later use we denote the value of x′ at the saddle point as x′‡. With x′ as the new reaction coordinate, we have

(3.66)

Inserting Eqs. (3.63) and (3.66) into Eq. (3.31), we indeed arrive at the Grote–Hynes result, Eq. (3.57c). Note that this rate constant can be written in the form

(3.67)

The above derivation amounts to choosing x′=x′‡ as the dividing surface of the multi-dimensional system of reaction coordinate plus thermal bath. Had we chosen x=x‡, based on the original reaction coordinate, as the dividing surface, we would have derived the standard transition-state theory result, Eq. (3.16). Pollak's use of a dividing surface based on the modified reaction coordinate can thus be seen as a form of variational transition-state theory. Berezhkovskii et al. (Reference Berezhkovskii, Pollak and Zitserman1992) have, in fact, carried out variational transition-state theory calculations using a planar dividing surface. They found that, when the harmonic approximation of Eq. (3.45) is invoked for the potential function, the optimized planar dividing surface is x′=x′‡. However, when anharmonicity is accounted for, the optimized planar dividing surface differs from x′=x′‡ in direction and no longer goes through the saddle point if U(x) is asymmetric with respect to x=x‡.

3.4 Solution to the turnover problem

As noted in subsection 3.2, a proper calculation of the rate constant at low frictions should account for a split boundary condition arising from the absorbing boundary around the bottom of the product well, which traps molecules moving to the left but not those moving to the right. Melnikov & Meshkov (Reference Melnikov and Meshkov1986) dealt with this split boundary condition by the Wiener–Hopf method. When the rate constant is written in the form

(3.68)

where k _+TST is given by Eq. (3.16), their result for the transmission coefficient is

(3.69)

In the last expression I _A‡ is given by Eq. (3.48), I _B‡ is the counterpart for the product well, and the function z(x) is given by

(3.70)

which goes to x when x→0 and to 1 when x→∞.

In the low-friction limit, Eq. (3.69) gives =γI _A‡I _B‡/(I _A‡+I _B‡)k _BT. This differs from the corresponding result of Kramers, γI _A‡/k _BT [see Eq. (3.47b)] by a factor I _B‡/(I _A‡+I _B‡), which arises from accounting for the fact that the orbit with total energy ΔU‡ traverses both the reactant and the product wells. In the high-friction limit, Eq. (3.69) predicts =1, i.e., the transition-state theory result. The validity of Eq. (3.69) is thus restricted to low to moderate frictions. Melnikov and Meshkov recognized that, when the transmission coefficient of Eq. (3.69) is multiplied to Eq. (3.54), which is Kramers' result for the rate constant at moderate to high frictions, the resulting expression

(3.71)

is numerically accurate for both the low to moderate and the moderate to high ranges of friction. This expression can thus be used for the full range of friction. Computer simulations have shown that Eq. (3.71) is indeed highly accurate for all frictions (Zhou, Reference Zhou1989) (see Fig. 3b).

Pollak et al. (Reference Pollak, Grabert and Hanggi1989) subsequently developed an improved solution to the turnover problem. Following the work of Pollak (Reference Pollak1986) outlined above, they modeled the thermal bath as a set of harmonic oscillators. The rate constant at low frictions was derived by studying the sampling of energy levels associated with the modified reaction coordinate x′, not those associated with x as done by Kramers and by Melnikov and Meshkov. When the friction is taken to the high range, this rate constant becomes identical to the result of Pollak (Reference Pollak1986), which was specifically derived for that range of friction. Therefore, the theory of Pollak et al. works for the full range of a friction with memory. Its predictions are found to agree well with results from computer simulations of Straub et al. (Reference Straub, Borkovec and Berne1986).

3.5 Multi-dimensional potential energy functions

A description restricted to a single reaction coordinate has limitations. Often it makes sense to explicitly consider other closely coupled coordinates. Below we present two perspectives on how to deal with motions on multi-dimensional energy surfaces in modeling rate processes.

3.5.1 Langer's generalization of Kramers' result and further developments

Following Kramers' work on a one-dimensional system, Langer (Reference Langer1969) studied the rate of escape from a multi-dimensional energy well through a saddle point. Consider first the case where motions along the coordinates, x_i, collectively denoted as x, are diffusive. Then the probability density ρ(x, t) satisfies the multi-dimensional Smoluchowski equation

(3.72)

where ρ_eq(x)∝exp[−U(x)/k _BT], and the diffusion matrix D is assumed to be position independent. [It should be noted that, for a polyatomic reactant molecule, the diffusion matrix can be very different in different parts of the conformational space (see, e.g., McCammon & Harvey, Reference Mccammon and Harvey1987). This positional dependence of D can significantly affect the reaction paths and rate constant, as already emphasized for diffusion along a one-dimensional reaction coordinate (subsection 3.2). In higher dimensions, the minimum-resistance path from the reactant well to the product well has been calculated for a position-dependent diffusion matrix by a variational approach (Berkowitz et al. Reference Berkowitz, Morgan, Mccammon and Northrup1983).] Around the bottom of the reactant well, located at x=x_A, a harmonic approximation is made:

(3.73a)

All the eigenvalues of the force-constant matrix e_A are assumed to be positive. Around the saddle point x=x‡, another harmonic approximation is made:

(3.73b)

One of the eigenvalues of the force-constant matrix e‡ is assumed to be negative and all others are positive.

Adopting Kramers' trick, Langer sought for a set of parallel planes near the saddle point, over each of which the splitting probability ζ_A(x)=ρ_ss(x)/ρ_eq(x) is a constant. This condition is satisfied when the normal of the planes is the eigenvector, s, corresponding to the negative eigenvalue, −ω₊‡, of the matrix (k _BT)⁻¹e‡⋅D. This negative eigenvalue is determined by

(3.74)

and the parallel planes are specified by the equation

(3.75)

The splitting probability around the saddle point is given by

(3.76)

where the magnitude of s is chosen so that

(3.77)

According to Eq. (3.76), the splitting probability has value 1/2 on the plane s⋅(x−x‡)=0. This plane, which goes through the saddle point, is thus the stochastic separatrix (see subsubsection 3.6.2 below). Note that, when the diffusion matrix is highly anisotropic, the eigenvector s, which is normal to the stochastic separatrix, may be very different from the eigenvector of the force-constant matrix e‡ corresponding to its negative eigenvalue.

Finally, the flux-over-population formulation gives the rate constant as

(3.78)

which has the same structure as Eq. (3.67). Specializing to a one-dimensional model, we recover Kramers' result [Eq. (3.46)] for high frictions. Interestingly, Berezhkovskii & Szabo (Reference Berezhkovskii and Szabo2005) showed that, if a one-dimensional reaction coordinate, x, is chosen along the eigenvector s in the multi-dimensional space, then Langer's result is identical to Kramers' one-dimensional result if the one-dimensional potential is calculated as the potential of mean force in x and the diffusion constant along x is chosen as s·D·s.

Now consider the case where motions along the coordinates are governed by the Langevin equation, which here takes the form

(3.79)

where m_i are the masses corresponding to coordinates x_i. Let the total number of coordinates be n. (For a polyatomic reactant molecule, overall translation and rotation correspond to zero eigenvalues of the force-constant matrix; such coordinates are assumed to be appropriately eliminated.) To simply notation, we absorb each mass into the corresponding coordinate and velocity; i.e., m_i ^1/2x_i and m_i ^1/2v_i will hereafter be written as x_i and v_i. The probability density ρ(x, v, t) satisfies the n-dimensional Fokker–Planck equation

(3.80)

where ρ_eq(x, v)∝exp[–(U(x)+(1/2)v⋅v)/k _BT].

Again, we seek for a set of parallel planes over which ζ(x, v)=ρ_ss(x, v)/ρ_eq(x, v) is uniform. [The splitting probability ζ_A(x, v) is ζ(x, −v).] The normal of these planes is the eigenvector corresponding to the negative eigenvalue, −ω₊‡, of the 2n×2n matrix

where 0 and I are n×n zero and identity matrices, respectively; f‡ is the n×n mass-scaled force-constant matrix, related to the original force-constant matrix e‡ via

(3.81)

and γ is the n×n matrix with elements γ_ij. The eigenvalue is determined by

(3.82)

Note that, at high frictions, the term is negligible, and Eq. (3.82) reduces to Eq. (3.74). Let the first and last n components of the eigenvector be denoted as s and t, respectively. The equation of the parallel planes is then

(3.83)

Around the saddle point, ζ(x, v) is given by

(3.84)

where the magnitudes of s and t are chosen so that

(3.85)

Finally, the flux-over-population formulation gives the rate constant as

(3.86)

where f_A is the n×n mass-scaled force-constant matrix at the bottom of the reactant well. Specializing to a one-dimensional model, we recover Kramers' result [Eq. (3.54)] for moderate to high frictions.

Berezhkovskii et al. (Reference Berezhkovskii, Pollak and Zitserman1992) considered the case where each coordinate of the multi-dimensional reactant system is coupled to a thermal bath consisting of a set of harmonic oscillators. Motion along each coordinate is then governed by a generalized Langevin equation [Eq. (3.79) with friction coefficients replaced by memory kernels]. They carried out variational transition-state theory calculations on the high-dimensional system of reactant plus thermal baths, using a planar dividing surface. Within the harmonic approximations of Eqs. (3.73), the rate constant is given by Eq. (3.86), but with the governing equation for ω₊‡ generalized to

(3.87)

where is a matrix with elements given by the Laplace transform of the memory kernels. A more cumbersome way to derive Eqs. (3.86) and (3.87), when motions on the multi-dimensional energy function are governed by generalized Langevin equations, is to first transform these coupled equations into an equation of motion for a single reaction coordinate, with a more complicated memory kernel, and then apply the Grote–Hynes theory (subsection 3.3) (Nitzan, Reference Nitzan1987).

As an aside, we note that the expression of the rate constant in terms of a time integral, Eq. (3.30), can be generalized to the present case of a multi-dimensional energy function. The result is (Bicout & Szabo, Reference Bicout and Szabo1997; Zhou, Reference Zhou2005a)

(3.88)

where θ_A(x) has value 1 if x is in the reactant well and 0 otherwise, and ρ_N(x, t |x₀) is the conditional probability density satisfying natural boundary conditions. This result will find use in the next section.

3.5.2 Agmon–Hopfield model

Agmon & Hopfield (Reference Agmon and Hopfield1983) considered the situation where an intramolecular reaction, such as electron transfer or geminate binding, is coupled to conformational fluctuations of the protein molecule. To model the situation, the approach presented so far would suggest that one treats explicitly the motion along the reaction coordinate x (e.g., the distance between the germinate ligand and its binding site) and accounts implicitly for the conformational fluctuations through their effects on the potential energy and the friction coefficient for x (see, e.g., D'Abramo et al. Reference D'Abramo, Di Nola and Amadei2009; Schaad et al. Reference Schaad, Zhou, Szabo, Eaton and Henry1993). Agmon and Hopfield took the opposite approach. The intramolecular reaction is implicitly modeled by a rate constant, k ₊, which depends on the protein conformation, here denoted with coordinate y. The conformational fluctuations are modeled explicitly as diffusion on an effective potential surface, U(y). On account of the reaction, the probability density ρ(y, t) satisfies the diffusion-reaction equation

(3.89)

where D is the effective diffusion constant. For a given normalized initial distribution ρ(y, 0), the quantity of interest is the survival probability

(3.90a)

Integrating both sides of Eq. (3.89) over y, we have

(3.90b)

In general, the rate coefficient, k(t), defined above is time dependent.

In the limit that diffusion along y is infinitely slow, the reactions at different y values proceed independently. Hence, as D→0,

(3.91a)

and the survival probability is

(3.91b)

In the opposite limit D→∞, equilibration among different y values occurs quickly, and ρ(y, t) is proportional to the equilibrium distribution, ρ_eq(y)=exp[−U(y)/k _BT]/∫dy exp[−U(y)/k _BT], but with a magnitude that decreases with increasing time. Then the rate coefficient k(t) in Eq. (3.90b) is simply the equilibrium average of k ₊(y):

(3.92a)

and

(3.92b)

becomes exponential. For intermediate values of D, Eq. (3.89) can only be solved numerically.

The simple model introduced in subsection 2.1 has the flavor of the Agmon–Hopfield model. Here motion (i.e., jumps) along the y-coordinate (i.e., discrete microstates) models the equilibration within the reactant state. Like the D→∞ limit just presented, fast intra-state equilibration leads to an exponential decay of the reactant population, and the decay constant is given by the equilibrium average of rate constants among the microstates [see Eq. (2.13b)].

3.6 Two useful quantities

3.6.1 Mean first passage time

As alluded to in section 2, the rate constant k ₊ is essentially the inverse of the mean first passage time from the reactant well to the product well. Formally, the mean first passage time, τ_FP, is the average time for a molecule to reach a boundary for the first time. That boundary is absorbing since the molecule never returns from it. τ_FP can be calculated from the survival probability, S(t), i.e., the probability that the molecule has not reached the absorbing boundary after a waiting time t. It is easy to see that the probability density, φ(t), of the first passage time t is . Hence,

(3.93)

Let the molecule be specified by coordinates collectively denoted as w (which could be just spatial coordinates x or both x and the corresponding velocities v). Suppose that the probability density ρ(w, t) is governed by

(3.94)

where (w) is an operator like those in Eqs. (3.72) and (3.80). The equilibrium distribution, ρ_eq(w), is a stationary solution. The survival probability S(t | w₀) starting from the initial position w₀ satisfies (Pontryagin et al. Reference Pontryagin, Andronov and Vitt1933)

(3.95)

The adjoint operator is given by (Risken, Reference Risken1989)

(3.96)

where Ɛw₀ means that the sign of each velocity component in w₀ is negated. Integrating both sides of Eq. (3.95) over time leads to

(3.97)

For the problem of diffusive motion on the one-dimensional energy surface of Fig. 3a, with an absorbing boundary at x=x _B, the solution of Eq. (3.97) is (Szabo et al. Reference Szabo, Schulten and Schulten1980)

(3.98)

Compared with Eq. (3.44), τ_FP(x _A) is almost identical to the inverse of Kramers' result for the rate constant, except that the upper limit of the second integral here is x, not x‡. Numerically, this difference hardly matters since the first integral is dominated by a region around x=x‡.

An elegant application of mean first passage time is in modeling diffusive motion on a rough potential surface (Zwanzig, Reference Zwanzig1988). Consider a potential function given by

where U ₀(x) is smooth and U ₁(x) is a rapidly oscillating perturbation. Over the range, Δx, in which U ₁(x) oscillates, U ₀(x) can be considered as a constant; hence,

Over a range much longer than Δx, the above integral can then be evaluated as the integral of . Hence,

Similarly,

with ψ⁺(x) defined analogous to ψ⁻(x). The mean first passage time of Eq. (3.98) now becomes

which is the same as if the molecule moves on an effective potential energy U ₀(x)+ψ^–(x) with an effective diffusion constant

(3.99)

Because , roughness in the potential always leads to a reduction in diffusion constant. This provides a simple explanation for the decrease in diffusion constant when a protein moves from three-dimensional diffusion in the bulk solution to one-dimensional diffusion along a DNA, as speculated previously (Berg et al. Reference Berg, Winter and Von Hippel1981; Schranner & Richter, Reference Schranner and Richter1978) [see subsubsection 4.3.7; an alternative, or additional proposed mechanism is that, while sliding along a groove of the DNA, the protein also rotates around the DNA axis, and this coupled rotation encounters significant solvent friction (Schurr, Reference Schurr1979)].

When U ₁ is random, and follows a Gaussian distribution with zero mean and variance ε, one finds ψ⁺(x)=ψ⁻(x)=ε²/2k _BT. Then the effective potential is changed by a trivial constant ε²/2k _BT, but the effective diffusion constant is reduced to Dexp[−(ε/k _BT)²], with the exponent having a quadratic dependence on temperature. Such a T ² dependence of the effective diffusion constant has been found for a random energy model of protein folding (Bryngelson & Wolynes, Reference Bryngelson and Wolynes1989).

3.6.2 Splitting probability

In subsection 3.2 and subsubsection 3.5.1, we introduced the splitting probability via its relation to a steady-state probability density. This relation can be generally written as

(3.100)

The steady-state probability density satisfies

(3.101)

with the boundary conditions

(3.102a)

(3.102b)

where x_A and x_B represent the bottoms of the reactant and product wells, respectively. It appears that the relation between ζ_A(w) and ρ_ss(w) has not been recognized previously. We now outline its proof.

The splitting probability ζ_A(w) satisfies the equation (Gardiner, Reference Gardiner1985)

(3.103)

with the boundary conditions

(3.104a)

(3.104b)

[If the roles of states A and B were reversed, Eqs. (3.103) and (3.104) would define the splitting probability, ζ_B(w), for the molecule to first reach the product state instead of the reactant state. It is clear that ζ_B(w)=1−ζ_A(w).] Combining Eqs. (3.96) and (3.101), it can be seen that ρ_ss(Ɛw)/ρ_eq(w) satisfies Eq. (3.103). In addition, from Eqs. (3.102), it can be seen that ρ_ss(Ɛw)/ρ_eq(w) satisfies the boundary conditions of Eqs. (3.104). Therefore, ρ_ss(Ɛw)/ρ_eq(w) and ζ_A(w) must be identical.

For the problem of diffusive motion on the one-dimensional energy surface of Fig. 3a, the splitting probability is

(3.105)

The integrals are dominated by a small region around x=x‡. Hence, ζ_A(x) is close to 1 from x=x _A up to a small distance to the left of x=x‡; it then undergoes a rapid transition to ~1/2 at x=x‡ and to 0 a small distance thereafter. Correspondingly, ρ_ss(x) undergoes a transition from ρ_eq(x) in the reactant well to 0 in the product well.

When the motion along x is governed by the Langevin equation, ρ_ss(x, v) undergoes a similar transition, along the new coordinate u=λ(x−x‡)−v. The value of ρ_ss(x, v) is ~1/2 at u=0, i.e., on the line v=λ(x−x‡). However, the splitting probability is ~1/2 on the line v=−λ(x−x‡). The sign switch on v is necessary because a molecule starting on the right of x=x‡ must have a leftward velocity for it to have the same splitting probability as a molecule starting at x=x‡ with a zero velocity.

In higher dimensions, the plane on which the splitting probability is 1/2 is known as the stochastic separatrix. Langer's solution to the rate problem shows that the stochastic separatrix is a plane that goes through the saddle point of the energy function; the orientation of the plane is determined by the force-constant matrix at the saddle point and the diffusion or friction matrix. In protein-folding kinetics, the splitting probability is given the name p _fold, and the stochastic separatrix is assumed to define the transition-state ensemble (Du et al. Reference Du, Pande, Grosberg, Tanaka and Shakhnovich1998).

3.7 Jump dynamics among discrete microstates

We now further consider the case where motion in the conformational space of a reactant molecule is represented by jumps among a set of discrete microstates. The occupational probabilities, ρ_i(t), are governed by a master equation:

(3.106)

where k_ji is the jump rate from microstate i to microstate j. The equilibrium distribution, ρ_{i eq}, is determined by the configurational integrals of the microstates. The jump rates satisfy the detailed balance condition:

(3.107)

In particular, master equations have been used to model protein-folding kinetics. In earlier work, the jump rates were made up, only constrained by the detailed balance condition (Alm et al. Reference Alm, Morozov, Kortemme and Baker2002; Cieplak et al. Reference Cieplak, Henkel, Karbowski and Banavar1998; Merlo et al. Reference Merlo, Dill and Weikl2005; Muñoz et al. Reference Muñoz, Henry, Hofrichter and Eaton1998; Schonbrun & Dill, Reference Schonbrun and Dill2003; Weikl et al. Reference Weikl, Palassini and Dill2004; Zwanzig, Reference Zwanzig1995); in more recent work (Buchete & Hummer, Reference Buchete and Hummer2008; Hummer & Kevrekidis, Reference Hummer and Kevrekidis2003), the jump rates are calculated from relatively short simulations of transitions between neighboring microstates [a similar idea is behind the so-called Markov state models (Chodera et al. Reference Chodera, Swope, Pitera and Dill2006, Reference Chodera, Singhal, Pande, Dill and Swope2007; Noe & Fischer, Reference Noe and Fischer2008)]. In matrix–vector form, Eq. (3.106) is

(3.108)

The jump matrix k has one zero eigenvalue, with ρ_eq as the corresponding eigenvector. The other eigenvalues, λ_l, l=1, 2, …, are positive; let λ_l be ordered from small to large.

Now suppose that the microstates are grouped into two states, A and B. When λ₁≪λ_l for all l>1, the transitions between the two states can be modeled well as rate processes, and the rate constants, k _±, are given by

(3.109)

The equilibrium occupational probabilities of the two states are

(3.110)

Instead of solving the eigenvalue problem, one can make a transition-state theory-type estimate for k _±; the presentation below follows the work of Zhou (Reference Zhou2008). As indicated by the derivation in subsubsection 3.1.1, a transition-state theory estimates k ₊ by the normalized total reaction flux from state A to state B, assuming that the occupation of the microstates in state A is according to the equilibrium distribution. That is,

(3.111)

Equation (2.13b) presented in subsection 2.1 can be viewed as a special case of Eq. (3.111). As the discussion of subsubsection 3.5.2 makes it clear, Eq. (3.111) becomes exact in the limit that intra-state equilibration is much faster than inter-state transition.

Like in the continuous case, k _+TST provides an upper bound for the rate constant, as has been shown explicitly. In the spirit of variational transition-state theory, one may therefore try to find the best k _+TST by varying the dividing surface between the reactant and product states. Here the dividing surface consists of all the allowed jumps (those with non-zero jump rates) between the reactant and product microstates. The optimal dividing surface is the one that minimizes the normalized reaction flux [Eq. (3.111)] and is hence referred to as the minimum-reaction-flux surface. [Note that, as recognized by Krivov & Karplus (Reference Krivov and Karplus2002), locating this surface is isomorphic to the minimum-cut problem in graph theory (Stoer & Wagner, Reference Stoer and Wagner1997).] Test against numerical solutions of master-equation models of protein folding shows that that the minimum-reaction-flux results work well. The minimum-reaction-flux surface separates microstates with large differences in p _fold and can be considered as constituting the transition-state ensemble in the discrete case. Among the forward jumps contributing to the minimum reaction flux, one can further identify the one(s) making the largest contribution as the ‘saddle point.’ However, the minimum-reaction-flux approach differs from traditional transition-state theory in one important respect: the dividing surface minimizing the normalized reaction flux is determined not only by the free energies of the microstates (as reflected by ρ_{i eq}) but also by the jump rates k_ji. In principle, dynamic information incorporated by the jump rates can lead to a dominant transition route (or minimum-resistance path) very different from that produced by the free-energy landscape alone. This situation is reminiscent of the continuous case with a position-dependent diffusion coefficient (subsubsection 3.5.1).

4.1 Formulation based on dissociation

One approach to calculate the association rate constant k _a is to first obtain the rate constant k _d for dissociation, which, as noted in subsection 2.2, is a unimolecular reaction, and then find k _a as the product of k _d and the association constant K _a [Eq. (2.17)] (Zhou, Reference Zhou2005a). This requires modeling the motion of the reactant pair in the native-complex energy well. To support results later in this section, here we model motion in the inner well and over the energy barrier as diffusive. We arrive at four useful results. The first is a clarification of the absorbing boundary condition usually used in modeling the diffusion-controlled regime. The second is a demonstration that, in a complete theory, the choice for the location of the transient complex does not affect the overall association rate constant. The third is a clarification of the reactivity that appears in the radiation boundary condition to account for a finite rate constant (i.e., k ₊) for barrier crossing into the native-complex energy well. The last is an expression for k _a based on a time integral.

The probability density in r, ρ(r, t), is assumed to satisfy the three-dimensional Smoluchowski equation

(4.5)

for the full range of the intermolecular distance r; here D is the diffusion constant for relative translation. Because of the centrosymmetry, the problem is effectively one dimensional, i.e., along r; the effect of the other two degrees of freedom is captured by the geometric factor 4πr ². We apply Kramers' result, Eq. (3.44), to obtain the dissociation rate constant. Setting the limits of the two integrals requires some care. One integral is over the range of the reactant state, which in the present case is the bound state; the lower limit is obviously r=0, and we choose the upper limit to be the top of the energy barrier, i.e., r=r‡. The other integral is from the bottom of the reactant well to the bottom of the product well; the lower limit is obviously r=a, and we choose the upper limit to be at infinite distance, since only then the reactant molecules are fully dissociated. The dissociation rate constant is thus

(4.6)

where we have appropriately inserted the geometric factor 4πr ².

The association constant for the present model is (Shoup & Szabo, Reference Shoup and Szabo1982; Zhou & Gilson, Reference Zhou and Gilson2009)

(4.7)

Therefore, the association rate constant is

(4.8)

This result is identical to what Debye (Reference Debye1942) found by imposing an absorbing boundary condition at r=a (see subsubsection 4.3.1). In the above derivation, there is no absorbing boundary condition at r=a; rather, a deep well is there. Therefore, the absorbing boundary condition used in modeling the diffusion-controlled regime corresponds to a deep and narrow energy well in front of a reflecting wall. This is the first main result of this subsection.

Note that the location of the transient complex does not appear in the expression for k _a. To make connection with the overall association rate constant given by Eq. (4.4), we write Eq. (4.8) as

The numerator of the first term can be identified with the inverse of the rate constant k ₊ of unimolecular transition to the native complex from the transient complex. The transient complex here is confined to the narrow region between r‡ and r‡+Δ (>R), with the latter serving as a reflecting boundary for the transient complex. The denominator of the first term can be identified with the association constant, , for the transient complex [cf. Eq. (4.7)]. The second term can be recognized as the diffusion-controlled rate constant, obtained by imposing an absorbing boundary condition at r=R:

(4.9)

Thus,

(4.10a)

This shows that k _a given by Eq. (4.8) can be written in the form of Eq. (4.4). Had we started with Eq. (4.4) and calculated the diffusion-controlled rate constant k _D and the activation-controlled rate constant k _A, we would have ended up with Eq. (4.8) as the final result for k _a, which, as noted above, is independent of the location of the transient complex. This demonstrates that, in a complete theory, the choice for the location of the transient complex does not affect the overall association rate constant.

Equation (4.10a) for k _a can be further written in a form obtained by imposing the radiation boundary condition of Collins & Kimball (Reference Collins and Kimball1949) at r=R (Shoup & Szabo, Reference Shoup and Szabo1982):

(4.10b)

Compared with the Collins–Kimball result [Eq. (4.18)], we find that the reactivity κ in the radiation boundary condition [Eq. (4.15a)] can be identified with k ₊Δ. Hence, κ is a way to implicitly model the transition from the transient complex to the native complex. This is the third main result of this subsection.

Finally, we apply the time-integral expression, Eq. (3.30), to find an alternative formula for the dissociation rate constant. In the present case, the ratio of the equilibrium occupation probabilities is

Note that the integral in the numerator essentially is the volume V of the container and is effectively infinite relative to the integral of the denominator. Hence ρ_Aeq→0 and ρ_Beq→1. Equation (3.30) specialized to the present case is thus

(4.11)

where ρ_N(r, t | r₀) is the conditional probability density satisfying natural boundary conditions. Combined with the equilibrium constant of Eq. (4.7), we arrive at the following expression for the association rate constant:

(4.12)

Using the conditional probability density to evaluate the above expression, we derive Eq. (4.8) as a good approximation for k _a when the inner well features a deep minimum. This expression of the rate constant is the last main result of this subsection, which will be further analyzed and extended below.

4.2 Activation-controlled regime: transition-state theory

While our focus is on the diffusion-controlled regime, for the sake of completeness, we mention the association rate constant predicted by the transition-state theory of Eyring (Reference Eyring1935). Corresponding to Eq. (3.31) for the rate constant of a unimolecular reaction, the rate constant for a bimolecular reaction is

(4.13)

where q _P and q _L are the partition functions of the protein and the ligand, respectively; is a partition function calculated for a protein–ligand pair without including the reaction coordinate (as signified by the prime), which is fixed at the value defining the transition state; U‡ is the protein–ligand interaction energy in the transition state; and V is the volume of the container.

To illustrate, consider the centrosymmetric model. We calculate the partition functions according to classical statistical mechanics. The protein molecule has only three degrees of freedom, for overall translation. Hence, q _P=h ⁻³V(2πk _BTm _P)^3/2, where m _P is the mass of the protein. Similarly, q _L=h ⁻³V(2πk _BTm _L)^3/2. involves five degrees of freedom: three for overall translation and two for relative translation orthogonal to r; the masses corresponding to these two types of motions are the total mass, m _P+m _L, and reduced mass, m _Pm _L/(m _P+m _L), respectively. Hence,

The activation-controlled rate constant is finally given by

(4.14)

As to be expected from a classical treatment, Planck's constant disappears from the final result.

4.3 Diffusion-controlled regime

We now explicitly model the overall translational and rotation diffusion of the reactant molecules and implicitly model the barrier crossing into the inner energy well. We show that the overall association rate constant can be written in the form of Eq. (4.4); thereafter, we focus on the diffusion-controlled rate constant k _D.

4.3.1 Centrosymmetric model

In the Smoluchowski (Reference Smoluchowski1917) approach to modeling diffusion-influenced bimolecular reactions, one works with the pair distribution function, P(r, t). It satisfies the Smoluchowski equation, Eq. (4.5), with a radiation boundary condition (Collins & Kimball, Reference Collins and Kimball1949)

(4.15a)

and a value

(4.15b)

The initial value of P(r, t) is

(4.15c)

corresponding to an equilibrium distribution. A time-dependent rate coefficient is defined as the total inward flux at r=R:

(4.16)

It can be shown that, if a protein binds irreversibly [as modeled by the radiation boundary condition of Eq. (4.15a)] with a ligand that is in excess, the decay of the protein concentration is governed by (Solc & Stockmayer, Reference Solc and Stockmayer1971; Szabo, Reference Szabo1989)

(4.17a)

or

(4.17b)

Note that C _p(t)/C_p(0) is equivalent to the survival probability of a single protein molecule surrounded by ligand molecules at concentration C _L.

The long-time limit, k _a(∞), of the rate coefficient introduced here corresponds to the association rate constant of subsection 4.1. We will continue to denote k _a(∞) simply as k _a. Solving for P(r,∞), one obtains

(4.18)

where k _D is the Debye result given by Eq. (4.9). This result for the overall association rate constant, as well as similar results later in this subsection, conforms to Eq. (4.4). As we explained above, the reactivity κ models the transition from the transient complex to the native complex. (When κ=0, the boundary becomes reflecting.) The limit κ→∞ corresponds to the situation where a reactant pair once reaching the transient complex is instantaneously transformed into the native complex. The rate of association is then limited by the relative translational diffusion to reach the transient complex; that is why the rate constant in this limit, k _D, is referred to as diffusion controlled. In this limit, the radiation boundary condition reduces to the absorbing boundary condition:

(4.19)

In short, the diffusion-controlled rate constant k _D is obtained by setting an absorbing condition at the transient complex. In the opposite limit κ→0, diffusion is able to replenish the reactant pairs lost by forming the native complex, so the initial equilibrium distribution is maintained. Then the steady-state value, k _a(∞), is the same as the initial value, k _a(0); i.e., the rate constant approaches

(4.20)

which is the activation-controlled limit.

When the interaction potential U(r) is absent, the Debye result [Eq. (4.9)] reduces to the well-known Smoluchowski formula

(4.21)

Throughout this section, we use a superscript ‘0’ to signify a rate constant calculated in the absence of an interaction potential; k _D⁰ is the limit of the rate constant set by random diffusion and is referred to as the basal rate constant. An attractive potential is expected to increase k _D. To illustrate, consider a Coulomb interaction potential, U(r)=−Q/r. Evaluating the Debye formula, we find

(4.22)

When the value of the potential at contact, U(R), is 9k _BT, the electrostatic attraction results in a modest increase of 9-fold in the association rate constant. This result will be contrasted below with much more significant electrostatic rate enhancement predicted from a more realistic model of protein–protein association.

The probability, S(t|r₀), that a protein–ligand pair, started at an initial separation r₀, has not formed the native complex at time t is known as the survival probability. It is related to the pair distribution function via

(4.23)

and hence satisfies the equation

(4.24)

with the boundary conditions

(4.25a)

(4.25b)

At t=∞, the left-hand side of Eq. (4.24) becomes 0, and the survival probability is called the escape probability. When the boundary at r=R is absorbing, the escape probability is very similar in form to the splitting probability of subsubsection 3.6.2. For the present model one has

(4.26)

which can be compared with the splitting probability given by Eq. (3.105). However, the behaviors of the splitting probability there and the escape probability here are very different. The energy function modeling the unimolecular transition between two states features a barrier, where the splitting probability undergoes a rapid transition from 0 to 1. In contrast, the energy function modeling intermolecular interaction has a longer range and is relatively monotonic; hence, the transition of the escape probability from 0 to 1 is more gradual. When no interaction potential is present, the escape probability is 1−R/r ₀.

As we clarified in subsection 4.1, an absorbing boundary is equivalent to a narrow deep energy well in front of a reflecting wall. If the absorbing boundary for calculating k _D is replaced by such an energy well, then the time-integral formula, Eq. (4.12), for the association rate constant can be used to calculate k _D. We thus have

(4.27)

where Δ, the width of the energy well, approaches zero, and the conditional probability density ρ_R(r, t|r₀) satisfies a reflecting boundary condition at r=R.Footnote ² From this formula for k _D, one can also derive the Debye result [Eq. (4.9)].

4.3.2 Anisotropic reactivity on protein molecule

As we emphasized, protein–ligand native complexes are stereospecific. The ligand binds to a specific site on the protein surface rather than the whole surface; hence, the reactant pair must satisfy constraints in both separation and relative orientation before the native complex can be formed. One then needs to model both relative translational diffusion and rotational diffusion of the reactant molecules. Here we treat the case where the ligand is point-like so that its orientation does not come into play; the case where rotational diffusion of the ligand is also modeled is deferred to the next subsubsection. For simplicity, we assume that the protein molecule is axially symmetric; we denote the axis of symmetric by the unit vector n_P. The pair distribution function now depends both on r and on n_P, and is denoted as P(r, n_P, t). On account of rotational diffusion, the governing equation for P(r, n_P, t) is

(4.28a)

where _P(n_P) is a rotational diffusion operator given by

(4.28b)

with D _P denoting the rotational diffusion constant of the protein molecule. The stereospecific binding site can be modeled by modifying the radiation boundary condition to one that is of the radiation type if the ligand falls inside the site and becomes reflecting otherwise (Solc & Stockmayer, Reference Solc and Stockmayer1971). That is, at r=R,

(4.29)

where is the unit vector along r, and θ* (, n_P) has value 1 if the angle between and n_P is within the angle span, δ_P, of the ‘reactive patch’ and value 0 otherwise.

The boundary condition of Eq. (4.29) is a mixed type (i.e., involving the pair distribution in one part of the boundary and its flux in another), which poses a mathematical challenge. To get around this problem, Shoup et al. (Reference Shoup, Lipari and Szabo1981) introduced an approximation. They assume that the flux is a constant over the reactive patch, and this constant is determined by requiring that Eq. (4.29a) is satisfied only on averaging over the reactive patch, not at every point over it. With this constant-flux approximation, the boundary condition at r=R now becomes

(4.30a)

where J (which depends on time for the time-dependent problem) is fixed by requiring

(4.30b)

with 〈· · ·〉* denoting an average over the reactive patch. In the absence of an interaction potential, the association rate constant is given by

(4.31a)

where

(4.31b)

(4.31c)

In the last two equations, F _P=(1−cosδ_P)/2 is the surface fraction covered by the reactive patch, μ_l=[l(l+1)D _P/D]^1/2R, P_l(x) are Legendre polynomials, and (π/2x)^1/2K_{l +1/2}(x) are modified Bessel functions of the third kind. The above result for k _a⁰ was essentially guessed by Solc & Stockmayer (Reference Solc and Stockmayer1973). Berg (Reference Berg1985) found that Λ_P/F _P can be approximated by

(4.31d)

where ν_P=[(1+D _PR ²/D)/2]^1/2. When D _P→0, Eq. (4.31c) reduces to

(4.31e)

The case of small reactive patches, corresponding to stereospecific binding, is of particular interest. A patch with angle span δ_P→0 on a spherical surface is equivalent to a disk with radius a=Rδ_P on an infinite reflecting plane. Applying the constant-flux approximation to the latter problem, Shoup et al. found that k _D⁰=(3π²/8)Da (assuming D _P=0). This is to be compared with the exact solution of the problem, k _D⁰=4Da (Hill, Reference Hill1975). These results show that the constant-flux approximation underestimates the diffusion-controlled rate constant by ~7% when δ_P→0. Note that, when δ_P→0, k _D⁰∝δ_P, not δ_P², as suggested by the surface fraction of the reactive patch. The weaker δ_p-dependence (and hence higher value) of k _D⁰ comes about due to the diffusive nature of the ligand molecule's motion. When a ligand molecule first reaches the reflecting part of the protein surface, it can come back and try again to reach the reactive patch.

Zhou (Reference Zhou1993) applied the constant-flux approximation to find the rate constant when a Coulomb potential U(r)=−Q/r is present. The result is given by

(4.32a)

where

(4.32b)

When the rotation diffusion constant D _P=0, the diffusion-controlled rate constant is

(4.32c)

where Q ₁=U(R)/2k _BT and (π/2x)^1/2I_{l +1/2}(x) are modified Bessel functions of the first kind. It was observed that, as δ_P→0, terms with larger and larger l become dominating in the sum over l; for those terms, the Q ₁-dependent factors become less and less important. Therefore, the whole denominator becomes nearly independent of Q ₁, and consequently

(4.33)

where k _D⁰ is the basal rate constant of Eq. (4.31e). The generality of Eq. (4.33) for small reactive patches (i.e., stereospecific binding) was recognized. In fact, this approximate formula was first suggested by results obtained by Brownian dynamics simulations, and the above deduction served to confirm the simulation results. Below we further elaborate on the generality of Eq. (4.33) and its use in modeling the effect of long-range electrostatic attraction on stereospecific protein–ligand binding.

4.3.3 Anisotropic reactivity on both protein and ligand molecules: free diffusion

Before we further consider the effect of the interaction potential, we present the more general binding model in which the protein molecule and the ligand molecule each bear a reactive patch, which for now is assumed to be axially symmetric. The pair distribution function, P(r, n_P, n_L, t), is now governed by

(4.34)

where _L(n_L) is the diffusion operator for ligand rotation [see Eq. (4.28b)]. The boundary condition at r=R becomes (Solc & Stockmayer, Reference Solc and Stockmayer1971)

(4.35)

where θ* (, n_P, n_L) has value 1 if the angles between and n_P and between and n_L are within the angle spans, δ_P and δ_L, of the respective reactive patches, and value 0 otherwise.

The rate constant for the one-patch model in the absence of an interaction potential is given by Eqs. (4.31a–c). To derive k _a⁰ for the present two-patch model, Solc & Stockmayer (Reference Solc and Stockmayer1973) first devised a reaction scheme for the one-patch model. The scheme involves two intermediates, in which the ligand is in contact with the reflecting part of the protein and with the reactive patch, respectively. The overall rate constant of this scheme was compared to Eqs. (4.31a–c) to identify the individual rate constants in the scheme. They then extended this scheme to the two-patch model, now involving four intermediates. The expression for k _a⁰ via this ‘quasi-chemical’ scheme for the two-patch model can again be written in the form of Eq. (4.31a), but now

(4.36a)

(4.36b)

where

(4.36c)

with Λ_P defined by Eq. (4.31c) and Λ_L defined analogously. For small patches, using Eq. (4.31d) for Λ_P and a similar approximation for Λ_L and keeping the lowest orders, one finds (Berg, Reference Berg1985)

(4.36d)

This result shows that, when δ_P and δ_L→0, k _D⁰∝δ_P²δ_L~δ_Pδ_L², not δ_P²δ_L², as suggested by the surface fractions of the reactive patches. Again, k _D⁰ is higher than expected from surface fractions because the diffusive nature of its motion allows the reactant pair to make repeated attempts to reach the absorbing boundary.

By applying the constant-flux approximation, Zhou (Reference Zhou1993) solved the diffusion equation for the two-patch model to find the time-dependent rate coefficient k _a⁰ (t). The result, in Laplace transform, is given by

(4.37a)

where k _a⁰ (0) is given by Eq. (4.36a). Note that the long-time limit of k _a⁰ (t) is given by the s→0 limit of s _a⁰ (s) and similarly for k _D⁰. The basal rate coefficient is

(4.37b)

where and is expressed in terms of Wigner's 3-j symbols:

The one-patch model is a special case with δ_L=π. In that case the only term in the summation over l ₂ is l ₂=0; it can be verified that the s→0 limit of Eq. (4.37b) reduces to Eq. (4.31c) of the one-patch model. Numerically, the results for the basal rate constant by the quasi-chemical scheme and by the constant-flux approximation are very close.

A third solution for the basal rate constant was obtained by Temkin & Yakobson (Reference Temkin and Yakobson1984), using a formula obtained by Doi (Reference Doi1975a, Reference Doib). Inspired by the closure approximation of Wilemski & Fixman (Reference Wilemski and Fixman1973), Doi derived a variational formula for the basal rate constant. With the simplest trial function, the basal rate constant for the present case takes the form

(4.38)

where, for notational simplicity, we have written (r, n_P, n_L) collectively as x, and ρ_R⁰ (x, t|x₀) is the conditional probability for free diffusion subject to a reflecting boundary condition at r=R. It can be recognized that Eq. (4.38) is a generalization of the time-integral formula of Eq. (4.27). Temkin and Yakobson's solution entailed finding the conditional probability and evaluating the integrals. The final result is identical to the s→0 limit of Eq. (4.37b), suggesting that the time-integral formula and the constant-flux approximation are equivalent. This equivalence was proven by Zhou & Szabo (Reference Zhou and Szabo1996a); see subsubsection 4.3.5.

The two-patch model was taken to its limit when the restriction on axial symmetry of the reactive patches was removed (Schlosshauer & Baker, Reference Schlosshauer and Baker2002). Now, for each reactant molecule, one additional rotational degree of freedom, i.e., rotation angle χ_P or χ_L around the unit vector n_P or n_L, needs to be accounted for. Using the constant-flux approximation, Schlosshauer and Baker were able to obtain the basal rate constant.

4.3.4 Anisotropic reactivity on both protein and ligand molecules: effect of interaction potential

Zhou (Reference Zhou1993) treated the two-patch model when a centrosymmetric interaction potential, U(r), is present. Applying the constant-flux approximation, he found

(4.39a)

The initial rate coefficient is

(4.39b)

The diffusion-controlled rate coefficient is

(4.39c)

where a prime signifies a derivative with respect to r, and f _l(r; μ) is a function that satisfies

(4.39d)

and decays to zero as r→∞. When U(r)=0, Eq. (4.39c) reduces to Eq. (4.37b).

While no solution of Eq. (4.39d) could be found for an arbitrary U(r), Zhou (Reference Zhou1997) was able to solve Eq. (4.39d) for two potential functions: a square-well potential and U(r)=−2k _BT ln[1+(exp(−U(R)/2k _BT)−1)R/r]. A main aim for finding explicit results of k _D for these potentials is to demonstrate the accuracy of the simple formula given by Eq. (4.33) for small patches. For example, at δ_P=δ_L=5°, the prediction of Eq. (4.33) overestimates the result given by Eq. (4.39c) for the logarithmic potential by merely 5% when U(R)=−9k _BT. The rate enhancement, , by this potential is ~8000. This dramatic rate enhancement is to be contrasted with the modest result of subsubsection 4.3.1 for the case where the reactant molecules have isotropic reactivity.

4.3.5 Generalization of Eq. (4.33)

Equation (4.33) presents an enormous simplification for accounting for the effect of interaction potentials on association rate constant. The individual results summarized above suggest that it is quite general for stereospecific binding. A general form of Eq. (4.33) was established through a pair of papers in 1996.

In the first paper, Zhou & Szabo (Reference Zhou and Szabo1996a) considered a general binding model in which the relative translational and overall rotational degrees of freedom of the reactant pair are collectively denoted as x. The transition from the transient complex to the native complex was modeled not by a radiation or absorbing boundary condition. Rather, the transient complex was assumed to be confined to a ‘reactive region’ in x space, specified by the function θ*(x) with value 1 in the reactive region and value 0 elsewhere. From the reactive region, the reactant pair forms the native complex with a rate constant k ₊. The governing equation for the pair distribution function is

(4.40)

where (x) is an operator modeling the diffusional motion of the reactant pair. This model is similar in spirit to the Agmon–Hopfield model [see Eq. (3.89)]; it can also be viewed as a microscopic implementation of reaction scheme (4.1). Because the transition from the transient complex to the native complex is accounted for by the reaction term in Eq. (4.40), the inner boundary is reflecting. The rate coefficient is given by

(4.41)

When the reactive region spans a patch on the inner boundary and a thickness Δ→0 in the normal direction, then the sink term in Eq. (4.40) is equivalent to a radiation boundary condition on the reactive patch, with κ given by k ₊Δ (which uncoincidentally is consistent with the discussion of subsection 4.1). The constant-flux approximation is derived when the pair distribution is assumed to follow the equilibrium distribution, P _eq(x)∝exp[−U(x)/k _BT], but with a time-dependent magnitude. Then it was shown that the rate coefficient can again be written as Eq. (4.39a), with

(4.42a)

(4.42b)

where 〈· · ·〉* denotes an average over the reactive region [the last step of Eq. (4.42a) applies to a small reactive region], and ρ_R (x, t | x₀) is the conditional probability satisfying

(4.43)

subject to a reflecting condition on the inner boundary. Equation (4.42b) reduces to the time-integral formulas for rate constants presented above in this subsection under different situations, including Eq. (4.38) for the two-patch model.

In a subsequent paper (Zhou, Reference Zhou1996), it was noted that, when the interaction potential U(x) is long ranged so that it is almost constant in the reactive region and its vicinity, then ρ_R (x, t | x₀) is almost independent of the potential, i.e., ρ_R (x, t | x₀)≈ρ_R⁰ (x, t | x₀), for x and x₀ both inside the reactive region. Consequently,

(4.44)

This is finally the general result that provides a powerful method for modeling the effect of long-range electrostatic attraction on stereospecific protein–ligand binding.

That a small reactive region is necessary for the validity of Eq. (4.44) can be demonstrated by the Debye result for the centrosymmetric model, in which the native complex has no requirement on relative orientation. In particular, Eq. (4.22), the Debye result for a Coulombic potential, shows a very different dependence on the potential than Eq. (4.44). As noted above, the rate enhancement given by Eq. (4.22) is much more modest than would be predicted by Eq. (4.44).

To complete the connection with reaction scheme (4.1), we note that the association constant for the transient complex is

(4.45a)

Combined with k _D given by Eq. (4.42b), we find that the rate constant for dissociating the transient complex by diffusion is

(4.45b)

The last result could also be derived from a time-integral formula, Eq. (3.88), of section 3. The right-hand side can be recognized as the mean residence time in the reactive region, which measures how long it takes for ligand molecules to move out of the reactive region via diffusion (while the pathway to form the native complex is turned off). The above derivation of Eq. (4.44) is the same as proving that the mean residence time is nearly independent of the long-range interaction potential. As observed by Zhou et al. (Reference Zhou, Wong and Vijayakumar1997), the smaller the reactive region, the faster the ligand population in it will decay, and hence the shorter the mean residence time. If the mean residence time is really short, then the ligand molecules that have diffused out will not move very far away. Now if the interaction potential does not vary significantly in the reactive region and its vicinity, then the ligand molecules would have moved in an essentially uniform potential. Consequently, k _−D would be nearly the same as if the potential is absent.Footnote ³

4.3.6 Atomistic models: transient-complex theory

The derivation of Eq. (4.44) makes it clear that its validity is subject to the condition that the reactive region is small while the interaction potential is long ranged. This condition is exactly met (fortuitously) for stereospecific protein–ligand binding under the influence of electrostatic attraction. The accuracy of Eq. (4.44) for atomistic protein–ligand models has been demonstrated by comparison against Brownian dynamics simulations (Zhou et al. Reference Zhou, Briggs and Mccammon1996, Reference Zhou, Wong and Vijayakumar1997, Reference Zhou, Briggs, Tara and Mccammon1998a).

The dependence of the diffusion-controlled rate constant on interaction potential given by Eq. (4.44) has an Arrhenius appearance. Such a dependence is expected for the activation-controlled rate constant k _A [see, e.g., Eq. (4.42a)], but totally unexpected for k _D [see, e.g., Eq. (4.22)]. Despite the Arrhenius appearance, k _D does not model an activation process.Footnote ⁴ The average interaction energy of the transient complex, 〈U〉*, which superficially resembles the activation energy in Eyring's and Kramers' theories, is typically negative (leading to rate enhancement) instead of being positive. Unlike the processes studied by Eyring and Kramers, where an energy barrier is what hinders the rate, the diffusion-controlled rate of protein–ligand binding is hindered by the translational and rotational constraints of the transient complex.

To calculate k _D, the transient complex has to be specified. Rather than being guided by any theoretical considerations, for many years the location of the transient complex (or, its equivalence in the form of an absorbing boundary) was proposed in an ad hoc way and often adjusted for best agreement with experiment (Altobelli & Subramaniam, Reference Altobelli and Subramaniam2000; Elcock et al. Reference Elcock, Gabdoulline, Wade and Mccammon1999; Gabdoulline & Wade, Reference Gabdoulline and Wade1997, Reference Gabdoulline and Wade2001; Gabdoulline et al. Reference Gabdoulline, Kummer, Olsen and Wade2003; Miyashita et al. Reference Miyashita, Onuchic and Okamura2004; Northrup & Erickson, Reference Northrup and Erickson1992; Northrup et al. Reference Northrup, Reynolds, Miller, Forrest and Boles1986; Spaar et al. Reference Spaar, Dammer, Gabdoulline, Wade and Helms2006; Vijayakumar et al. Reference Vijayakumar, Wong, Schreiber, Fersht, Szabo and Zhou1998; Zhou, Reference Zhou1993; Zhou et al. Reference Zhou, Briggs and Mccammon1996, Reference Zhou, Wong and Vijayakumar1997, Reference Zhou, Briggs, Tara and Mccammon1998a). In order to predict association rate constants from theory alone, the transient complex has to be specified without reference to experiment. The problem is especially challenging for protein–protein association, given the severity of orientational constraints on both reactant molecules. A solution was proposed by Alsallaq and Zhou (Alsallaq & Zhou, Reference Alsallaq and Zhou2007a), based on analyzing the interaction energy landscape of associating proteins.

As discussed in subsection 4.1, in a complete theory, the overall association rate constant k _a should not be sensitive to where the transient complex is placed. If it is placed far away from the native complex, then k _D will be large but k ₊ will be small. Conversely, if it is placed very close to the native complex, then k _D will be reduced but k ₊ will become very large. Either way, Eq. (4.4) is expected to give nearly the same result for k _a. However, given the considerable difficulty and uncertainty in the calculation of k ₊, it is highly desirable to use k _D as a close approximation for k _a. Then there is an optimal location for placing the transient complex (Zhou, Reference Zhou2001a). If it is placed too far from the native complex, then the resulting k _D would not be a useful approximation for k _a. On the other hand, placing the transient complex too close to the native complex would mean that short-range interactions and conformational rearrangement have to be dealt with in calculating k _D. The native complex sits in a deep well in the interaction energy landscape. The optimal placement for the transient-complex ensemble is at the outer boundary of the native-complex energy well (Alsallaq & Zhou, Reference Alsallaq and Zhou2007a; Zhou, Reference Zhou2001a) (Fig. 4b).

The algorithm for identifying the transient complex was based on the following observation: inside the native-complex energy well, translation and rotation are restricted, but, once outside, the proteins gain significant translational and rotational freedom (Zhou, Reference Zhou2001a) (Fig. 5a). Thus, the outer boundary of the native-complex energy well coincides with the onset of translational and rotational freedom. To simplify the calculations required for determining the transient complex, the short-range interaction energy stabilizing the native complex was modeled by the number of contacts, N _c, formed between the protein partners. Translational and rotational freedom was measured by σ_χ(N _c), the standard deviation of the rotation angle χ in configurations with a given contact level N _c. A sharp increase in σ_χ with decreasing N _c marks the onset of translational and rotational freedom and hence the location of the transition complex (Fig. 5b).

Fig. 5.

Specification of the transient complex for the barnase–barstar protein pair (Alsallaq & Zhou, Reference Alsallaq and Zhou2008). (a) Scatter plot of allowed (i.e., clash-free) configurations. Each scatter point represents a cluster of allowed configurations with the indicated contact number (N _c) and angle (χ) of relative rotation. The N _c level defining the transient complex is shown in dark color. (b) Transition of the standard deviation of χ, σ_χ, from the native complex (with high contact numbers) to the unbound state (with low contact numbers). The start of the sharp increase in σ_χ, as indicated by an arrow, marks the transient complex.

Since k _D is used as the prediction for k _a, Eq. (4.44) can be rewritten as an equation for k _a:

(4.46)

where k _a⁰ now denotes the basal rate constant, i.e., the rate constant for reaching the transient complex by translational and rotational diffusion in the absence of any biasing force, and the interaction potential is denoted with a subscript ‘el’ to signify that only electrostatic interactions are to be included. The neglect of short-ranged non-electrostatic effects from the Boltzmann factor is based on two considerations. First, transient-complex configurations are separated by at least one layer of solvent (Alsallaq & Zhou, Reference Alsallaq and Zhou2007b, Reference Alsallaq and Zhou2008); therefore, short-ranged forces such as hydrophobic and van der Waals interactions are relatively weak in the diffusion process, leading to the transient complex. Second, short-range interactions, even when present within the transient complex, contribute much less to rate enhancement (i.e., k _a/k _a⁰) compared to long-range interactions (see next subsubsection). However, short-range interactions are essential for determining the location and size of the transient-complex ensemble in configurational space, which in turn affect the magnitude of k _a⁰. A transient-complex ensemble that is less restricted in translation and rotation will lead to a higher k _a⁰.

The algorithm for identifying the transient-complex ensemble presented above along with Eq. (4.46) constitutes the transient-complex theory. The theory provides an interesting dissection on the wide spectrum of observed rate constants of protein–protein association (Alsallaq & Zhou, Reference Alsallaq and Zhou2008). Brownian dynamics simulations have found that basal rate constants, set by unbiased diffusion, are in the range of 10⁴–10⁶ M⁻¹ s⁻¹ (Alsallaq & Zhou, Reference Alsallaq and Zhou2008; Qin & Zhou, Reference Qin and Zhou2009). This is the range of rate constants observed on antibody–antigen binding (Foote & Eisen, Reference Foote and Eisen1995; Hoffman et al. Reference Hoffman, Labranche, Zhang, Canziani, Robinson, Chaiken, Hoxie and Doms1999; Wassaf et al. Reference Wassaf, Kuang, Kopacz, Wu, Nguyen, Toews, Cosic, Jacques, Wiltshire, Lambert, Pazmany, Hogan, Ladner, Nixon and Sexton2006); presumably, antibodies are not specifically optimized for fast binding, and antibody–antigen binding processes are described well as unbiased diffusion. However, many proteins are observed to associate faster than the basal rate constant (Baerga-Ortiz et al. Reference Baerga-Ortiz, Rezaie and Komives2000; Candia et al. Reference Candia, Garcia and Latorre1992; Darling et al. Reference Darling, Kuchibhotla, Glaesner, Micanovic, Witcher and Beals2002; Escobar et al. Reference Escobar, Root and Mackinnon1993; Gianni et al. Reference Gianni, Engstrom, Larsson, Calosci, Malatesta, Eklund, Ngang, Travaglini-Allocatelli and Jemth2005; Goldstein & Miller, Reference Goldstein and Miller1993; Hemsath et al. Reference Hemsath, Dvorsky, Fiegen, Carlier and Ahmadian2005; Johnson et al. Reference Johnson, Mccoy, Bingman, Phillips and Raines2007; Korennykh et al. Reference Korennykh, Piccirilli and Correll2006; Miller, Reference Miller1990; Murrell-Lagnado & Aldrich, Reference Murrell-Lagnado and Aldrich1993; Park & Raines, Reference Park and Raines2001; Radic et al. Reference Radic, Kirchhoff, Quinn, Mccammon and Taylor1997; Schreiber & Fersht, Reference Schreiber and Fersht1996; Shapiro et al. Reference Shapiro, Ruiz-Gutierrez and Chen2000; Shen et al. Reference Shen, Hage and Sebald1996; Stewart & Van Bruggen, Reference Stewart and Van Bruggen2004; Terlau et al. Reference Terlau, Shon, Grilley, Stocker, Stuhmer and Baldomero1996; Uter et al. Reference Uter, Gruic-Sovulj and Perona2005; Walker et al. Reference Walker, Moore, James and Kleanthous2003; Wallis et al. Reference Wallis, Moore, James and Kleanthous1995; Wendt et al. Reference Wendt, Leder, Harma, Jelesarov, Baici and Bosshard1997). These association rate constants must have been enhanced by a biasing force. The electrostatic surfaces of the reactant partners in these systems indeed show complementarity (Alsallaq & Zhou, Reference Alsallaq and Zhou2008; Schreiber et al. Reference Schreiber, Haran and Zhou2009), indicating that electrostatic attraction provides the rate enhancement. The transient-complex theory suggests that the narrow range of 10⁴–10⁶ M⁻¹ s⁻¹ serves as a demarcation: lower rate constants implicate activation control, whereas rate constants in this range and higher implicate diffusion control, with the higher rate constants further implicating electrostatic attraction.

As recognized previously (Zhou, Reference Zhou2001a, Reference Zhou2003a), the transient-complex theory provides a nice explanation for a widely observed phenomenon regarding the effects of ionic strength on protein association kinetics. For many protein complexes that apparently are under diffusion control, the association and dissociation rate constants show disparate dependences on ionic strength: k _a decreases significantly with increasing ionic strength, whereas k _d is only modestly affected by ionic strength (Baerga-Ortiz et al. Reference Baerga-Ortiz, Rezaie and Komives2000; Candia et al. Reference Candia, Garcia and Latorre1992; Darling et al. Reference Darling, Kuchibhotla, Glaesner, Micanovic, Witcher and Beals2002; Escobar et al. Reference Escobar, Root and Mackinnon1993; Gianni et al. Reference Gianni, Engstrom, Larsson, Calosci, Malatesta, Eklund, Ngang, Travaglini-Allocatelli and Jemth2005; Goldstein & Miller, Reference Goldstein and Miller1993; Hemsath et al. Reference Hemsath, Dvorsky, Fiegen, Carlier and Ahmadian2005; Miller, Reference Miller1990; Murrell-Lagnado & Aldrich, Reference Murrell-Lagnado and Aldrich1993; Radic et al. Reference Radic, Kirchhoff, Quinn, Mccammon and Taylor1997; Schreiber & Fersht, Reference Schreiber and Fersht1993; Shen et al. Reference Shen, Hage and Sebald1996; Stewart & Van Bruggen, Reference Stewart and Van Bruggen2004; Walker et al. Reference Walker, Moore, James and Kleanthous2003; Wallis et al. Reference Wallis, Moore, James and Kleanthous1995; Wendt et al. Reference Wendt, Leder, Harma, Jelesarov, Baici and Bosshard1997). Ionic strength serves to modulate the magnitude of the electrostatic attraction between two proteins; the attraction is weakened when ionic strength is increased, leading to significant decrease in k _a. Under diffusion control, k _d≈(k ₋/k ₊)k _−D [Eq. (4.2b)]; neither k _−D nor k ₋/k ₊ has much dependence on ionic strength. The former is because k _−D does not depend much at all on the electrostatic attraction (see preceding subsubsection). The latter comes about because of the close proximity of the transient and the native complexes; then ionic strength modulates electrostatic interactions in them to similar extents and consequently has little effect on their equilibrium constant k ₋/k ₊.

Quantitative test against experimental data on many protein–protein and protein–RNA complexes (Law et al. Reference Law, Linde, Chambers, Oubridge, Katsamba, Nilsson, Haworth and Laird-Offringa2006; Radic et al. Reference Radic, Kirchhoff, Quinn, Mccammon and Taylor1997; Schreiber & Fersht, Reference Schreiber and Fersht1993, Reference Schreiber and Fersht1996; Shen et al. Reference Shen, Hage and Sebald1996; Wallis et al. Reference Wallis, Moore, James and Kleanthous1995) has clearly demonstrated the predictive power of the transient-complex theory (Alsallaq & Zhou, Reference Alsallaq and Zhou2007b, Reference Alsallaq and Zhou2008; Qin & Zhou, Reference Qin and Zhou2008). In a recent study, the theory was used to dissect a record-setting rate constant, at >10¹⁰ M⁻¹ s⁻¹ (Korennykh et al. Reference Korennykh, Piccirilli and Correll2006), of a ribotoxin binding to a biologically essential RNA loop on the ribosome (Qin & Zhou, Reference Qin and Zhou2009). Electrostatic attraction provides a 5×10⁶-fold rate enhancement. In comparison, the binding rate constant to the isolated RNA loop is electrostatically enhanced by 5×10²-fold, in line with the results found in other protein–protein and protein–RNA complexes (Alsallaq & Zhou, Reference Alsallaq and Zhou2007b, Reference Alsallaq and Zhou2008; Qin & Zhou, Reference Qin and Zhou2008). There are two contributions to the additional 10⁴-fold rate enhancement. First, the rest of the ribosome provides extra electrostatic attraction to the ribotoxin. Second, neighboring ribosomal proteins reshape the binding interface to position the transient complex into a region in configurational space where the electrostatic attraction between the ribotoxin and the RNA loop is particularly strong.

4.3.7 Facilitation by non-specific binding

A biological problem of great interest is the search by DNA-binding proteins for specific sites on genomic DNA. This potentially is a slow process, akin to finding a needle (i.e., a specific site) in a haystack (i.e., genomic DNA). That it is completed in short times routinely in cells prompted Adam & Delbruck (Reference Adam, Delbruck and Davidson1968) to look for an explanation. Their key observation is that, by non-specifically binding to the DNA surface, a protein searches for a specific site in a one-dimensional space rather than three-dimensional spaces. The search time, as their calculations of first mean passage time show, in the lower dimensionality is significantly shorter. This reduction of dimensionality idea has been further developed (Berg et al. Reference Berg, Winter and Von Hippel1981; Berg & Ehrenberg, Reference Berg and Ehrenberg1982; Halford & Marko, Reference Halford and Marko2004; Richter & Eigen, Reference Richter and Eigen1974; Riggs et al. Reference Riggs, Bourgeois and Cohn1970; Schranner & Richter, Reference Schranner and Richter1978). In particular, Berg et al. (Reference Berg, Winter and Von Hippel1981) phenomenologically described four modes of non-specific binding-facilitated translocation along the DNA: sliding, hopping, jumping, and intersegment transfer. Berg & Ehrenberg (Reference Berg and Ehrenberg1982) empirically incorporated diffusion on the DNA surface, or surface diffusion, into the Smoluchowski approach to diffusion-influenced reactions.

Zhou & Szabo (Reference Zhou and Szabo2004) introduced a more fundamental and realistic treatment, in which non-specific binding is accounted for by a short-range attractive potential around the DNA surface. The first three modes described by Berg et al. (Reference Berg, Winter and Von Hippel1981) are all encompassed in this treatment and there is no need (and rigorously it is impossible) to distinguish them. (The fourth mode, intersegment transfer, applies to bivalent proteins and is not dealt with in this work.) Zhou and Szabo modeled the protein as an isotropically reactive sphere, the DNA as an infinite cylinder (with contact distance R), and the specific binding site as a reactive strip (with height 2η); the short-range attractive potential, U(r), was assumed to be axially symmetric (r here denotes distance to the cylinder axis). Applying the constant-flux approximation, the rate coefficient again can be written in the form of Eq. (4.39a). The diffusion-controlled rate coefficient is

(4.47a)

where μ=(λ²+s/D)^1/2R, and f(r; μ) satisfies

(4.47b)

and decays to zero as r→∞.

For a square-well potential with width Δ and depth −U*, solution of Eq. (4.47b) leads to

(4.47c)

where

(4.47d)

with μ₁=μ(1+Δ/R). For a model with R=30 Å and η=3 Å, at U*=−9k _BT, Eqs. (4.47) predict rate enhancements of 18- and 38-fold, respectively, for Δ=1 and 5 Å. These enhancements are much more modest than if k _D scales with exp(−U*/k _BT) [Eq. (4.44)], due to the short-range nature of such an interaction potential. On the other hand, the non-specifically bound complex is probably also stabilized by non-electrostatic interactions, so the magnitude of the contact potential is much greater than the value from electrostatic attraction alone.

The above treatment was extended to account for the fact that the diffusion constant for one-dimensional diffusion along the DNA is much smaller than that for three-dimensional diffusion in the bulk solution (Blainey et al. Reference Blainey, Van Oijen, Banerjee, Verdine and Xie2006; Laurence et al. Reference Laurence, Kwon, Johnson, Hollars, O'Donnell, Camarero and Barsky2008; Wang et al. Reference Wang, Austin and Cox2006), and to the case of finite DNA lengths (Zhou, Reference Zhou2005b). We note that similar modeling has been used to study the role of non-specific binding in the search for membrane-bound receptors (Adam & Delbruck, Reference Adam, Delbruck and Davidson1968; Berg, Reference Berg1985; Richter & Eigen, Reference Richter and Eigen1974; Zhou & Szabo, Reference Zhou and Szabo2004).

4.4 Influence of conformational fluctuations

So far, we have assumed that the reactant molecules are frozen in their native conformations. How do conformational fluctuations affect the association rate constant? Szabo et al. (Reference Szabo, Shoup, Northrup and Mccammon1982) addressed this question by introducing a simple form of conformational fluctuations into the single-patch model discussed above. They assumed that the ligand molecule can alternate between two conformations: an ‘active’ conformation for which binding can proceed just as have been modeled through a radiation boundary condition on a reactive patch and an ‘inactive’ conformation for which the reactive patch is reflecting. The two conformations switch via rate processes:

(4.48)

The conformational switch, referred to as gating, is assumed to be uncoupled to the translational and rotational diffusion of the reactant molecules. Again using the constant-flux approximation, Szabo et al. were able to express the ‘gated’ association rate constant, k _G, to the Laplace transform of the ‘ungated’ rate coefficient:

(4.49)

where ω₀=ω₀₊+ω₀₋. This result has been extended to the time-dependent rate coefficient k _G(t) (in Laplace space) and to the case where both reactant molecules are gated (Zhou & Szabo, Reference Zhou and Szabo1996a).

When the gating is slow (i.e., ω₀₊ and ω₀₋→0), ω_0a(ω₀)→k _a; Eq. (4.49) becomes

(4.50a)

In this limit, a ligand molecule can contribute to the inward flux on the inner boundary only if it is initially in the active conformation, which it adopts with the equilibrium probability ρ_+eq. Hence, the total rate constant is this equilibrium probability times the rate constant for a ligand molecule that stays in the active form. In the opposite limit (i.e., ω₀₊ and ω₀₋→∞), ω_0a(ω₀)→k _a(0); therefore,

(4.50b)

where we have expressed k _a in terms of k _D and k _a(0) [see Eq. (4.32a)]. This is the rate constant for ungated binding with effective reactivity ρ_+eqκ; that effective reactivity is just what is expected when the conformational switch is fast. If the binding of the active form is diffusion controlled, i.e., κ→∞, Eq. (4.50b) reduces to k _G≈k _D. That is, the diffusion-controlled rate constant is unhindered when the conformational switch is fast. This result is very significant. It means that diffusion-controlled rate constants calculated under the assumption that the reactant molecules are frozen in their native conformations, as is the case for the transient-complex theory, are valid as long as the conformational fluctuations of the reactant molecules are fast.

The calculation of k _G(t), based on a single protein–ligand pair, is the same whether the protein or ligand is gated. However, the solution to the many-body problem of a protein molecule surrounded by an excess of ligand molecules does depend on which reactant is gated (Zhou & Szabo, Reference Zhou and Szabo1996b). When the ligand molecules are gated, effectively they independently bind with the protein molecule, and Eq. (4.17b) still holds for the survival probability of the protein, except that the rate coefficient is now replaced by k _G(t):

(4.51a)

However, when the protein molecule is gated, the dynamics of the ligand molecules around it become coupled: when the protein undergoes an active to inactive switch, all the surrounding ligand molecules simultaneously sense a change in the boundary condition from radiation to reflection. Zhou & Szabo (Reference Zhou and Szabo1996b) proposed the following rate equations:

(4.51b)

for the probabilities, _c(t), that the protein survives at time t while in conformation c, c=‘a’ for active and ‘i’ for inactive, respectively; (t)=_a(t)+_i(t). The solution of these rate equations is in good agreement with computer simulations. Under fast gating, Eq. (4.51b) reduces to Eq. (4.51a).

Some enzymes have buried active sites, with narrow tunnels leading to the exterior; conformational fluctuations of bottlenecks, or gates, along the tunnels must occur to allow for substrate access and product release (Zhou & McCammon, Reference Zhou and Mccammon2010). When a gate along the tunnel to the active site switched between closed and open states according to scheme (4.48), an approximate result for the gated binding rate constant similar to Eq. (4.49) was found (Zhou, Reference Zhou1998):

(4.52)

where k _a again is the ungated rate constant and J(t) is the total flux across the gate at time t when the substrate was started from an equilibrium distribution confined to the outside of the gate. Ĵ(s) has the following limiting values:

(4.53a)

(4.53b)

where A _g is the area of the gate and 〈· · ·〉_g denotes an average over the cross-section of the gate. Using these results, we can find the gated rate constant in the slow gating limit:

(4.54a)

and in the fast gating limit:

(4.54b)

The last result states that, under fast gating, the gated rate constant approaches the ungated rate constant, even if binding when the gate is open is not limited by diffusion.

Zhou et al. (Reference Zhou, Wlodek and Mccammon1998b) used Eq. (4.52) to study the substrate binding rate constant to acetylcholinesterase, an enzyme with a buried active site but still affording one of the highest k _cat/K _M values, at ~10⁹ M⁻¹ s⁻¹ (Nolte et al. Reference Nolte, Rosenberry and Neumann1980; Ordentlich et al. Reference Ordentlich, Barak, Kronman, Flashner, Leitner, Segall, Ariel, Cohen, Velan and Shafferman1993; Pryor et al. Reference Pryor, Selwood, Leu, Andracki, Lee, Rao, Rosenberry, Doctor, Silman and Quinn1992; Vellom et al. Reference Vellom, Radic, Li, Pickering, Camp and Taylor1993). Analysis of the gating dynamics obtained from molecular dynamics simulations showed that, for the intended substrate, acetylcholine, it is close to the fast-gating regime, and the gated rate constant is lower than the ungated rate constant by just 2-fold, even though ρ_+eq is only 2%. Combined with the ungated rate constant previously calculated from Brownian dynamics simulations (Zhou et al. Reference Zhou, Briggs and Mccammon1996), the gated rate constant is found to be in the range of measured k _cat/K _M values. However, for a somewhat larger substrate, modeling butyrylcholine, the gated rate constant is three orders of magnitude lower (ρ_+eq decreases to 8×10⁻⁵), in line with similar decreases in k _cat/K _M observed for that substrate (Ordentlich et al. Reference Ordentlich, Barak, Kronman, Flashner, Leitner, Segall, Ariel, Cohen, Velan and Shafferman1993; Pryor et al. Reference Pryor, Selwood, Leu, Andracki, Lee, Rao, Rosenberry, Doctor, Silman and Quinn1992; Vellom et al. Reference Vellom, Radic, Li, Pickering, Camp and Taylor1993). Conformational gating thus affords the enzyme exquisite substrate selectivity. That gating is essential in the selectivity is supported by the observation that when bulky side chains constituting the gate are mutated into smaller ones, presumably bringing the gating dynamics for butyrylcholine toward the fast-gating regime through an increase in ρ_+eq, k _cat/K _M for butyrylcholine is significantly increased (Ordentlich et al. Reference Ordentlich, Barak, Kronman, Flashner, Leitner, Segall, Ariel, Cohen, Velan and Shafferman1993; Pryor et al. Reference Pryor, Selwood, Leu, Andracki, Lee, Rao, Rosenberry, Doctor, Silman and Quinn1992; Vellom et al. Reference Vellom, Radic, Li, Pickering, Camp and Taylor1993).

4.5 Intermolecular versus intramolecular binding

When the protein and ligand are covalently linked, their binding becomes intramolecular. Under an approximate treatment of the linker, a simple relation between the equilibrium constants of the inter- and intramolecular binding processes has been derived (Zhou, Reference Zhou2001b; Zhou & Gilson, Reference Zhou and Gilson2009). The linker is assumed to be flexible and to have negligible interactions with the protein and the ligand; the end-to-end distance at which the probability density, ρ_L(r), of the end-to-end vector r is maximal is assumed to be much greater other relevant lengths, such as the size of the native-complex energy well, and the range of a long-range potential if present. Then

(4.55)

Throughout this subsection we use a superscript L to denote quantities for the intramolecular process; a denotes the end-to-end vector in the native complex. What are the corresponding relations between the rate constants of the inter- and intramolecular binding processes?

When the binding is activation controlled, it is easy to see that

(4.56a)

(4.56b)

When the binding is diffusion controlled, these relations for the rate constants still hold if K _a^L=K _aρ_L(a)≪1 and the linker does not affect the relative diffusion constant of the protein and ligand. The justification is the same as the one leading to Eq. (4.46), which captures the effect of long-range electrostatic attraction on k _a (Zhou, Reference Zhou2002, Reference Zhou2001c). A check of Eq. (4.56a) is provided by the rate constant for the diffusion-controlled end-to-end contact formation of a polymer chain obeying Gaussian statistics, obtained by Szabo et al. (Reference Szabo, Schulten and Schulten1980) by calculating the mean first passage time. Let the contact distance be a and the mean square of the end-to-end distance be 〈r ²〉. The probability density for the end-to-end vector is (3/2π〈r ²〉)^3/2exp(−3r ²/2〈r ²〉). When a≪〈r ²〉^1/2, which is equivalent to K _a^L≪1, Szabo et al.'s result is

The first factor, 4πDa, can be recognized as the diffusion-controlled rate constant for intermolecular contact formation; the second factor, (3/2π〈r ²〉)^3/2, can be recognized as the probability density of the end-to-end vector at r=a [when a≪〈r ²〉^1/2, exp(−3a ²/2〈r ²〉)→1]. Hence, Szabo et al.'s result conforms to Eq. (4.56a).

Equation (4.55) has been used to derive the association constant for a bivalent ligand, which features two fragments connected by a linker; the framents bind to two separate sites (named 1 and 2) on the protein (Zhou, Reference Zhou2001b). After binding the first fragment, the binding of the second fragment becomes intramolecular. The overall association constant of the bivalent ligand is thus

(4.57)

Here we use a superscript B to denote quantities for the bivalent ligand; K _a1 and K _a2 are the association constant for the first fragment when the second is absent and vice versa; and a again denotes the end-to-end vector of the linker in the native complex.

The native complex can be formed by first occupying site 1 and then site 2 or by proceeding in the reverse order:

Making a steady-state approximation for each singly bound species, the overall rate constant for binding the bivalent ligand is

(4.58a)

Applying the approximation of Eq. (4.56a) to k _a1^L and k _a2^L and expressing k _d1 and k _d2 as k _a1/K _a1 and k _a2/K _a2, respectively, we find (Zhou, Reference Zhou2003b)

(4.58b)

However, validity of Eq. (4.56a) requires that K _a1ρ(a) and K _a2ρ(a)≪1; hence,

(4.58c)

4.6 Reversibility

Up to now we have treated the protein–ligand binding process as though it is irreversible. The kinetics of reversible diffusion-influenced binding has long attracted the attention of theorists (Agmon, Reference Agmon1984; Gopich & Doktorov, Reference Gopich and Doktorov1996; Gopich & Szabo, Reference Gopich and Szabo2002; Lee & Karplus, Reference Lee and Karplus1987; Lukzen et al. Reference Lukzen, Doktorov and Burshtein1986; Naumann, Reference Naumann1994; Szabo, Reference Szabo1991; Yang et al. Reference Yang, Lee and Shin1998). For simplicity we only consider the pseudo-first-order limit, where effectively a single protein molecule is surrounded by excess ligand molecules. Upon binding with one ligand molecule, the protein surface becomes reflecting for all the other diffusing ligand molecules. The fact that the binding of a single ligand molecule instantaneously changes the boundary condition for the diffusion of all the other ligand molecules renders the problem many-body in nature. In general the kinetics is quite complicated, but under two limiting situations it is simplified. The first is when diffusion is infinitely fast (and the binding is thus activation controlled). Then a newly released ligand molecule will immediately equilibrate with the other ligand molecules to reach the equilibrium distribution, and consequently the ordinary rate equation, Eq. (2.15), applies. The relaxation function Y(t)=[C _P(t)−C _Peq]/[C _P(0)−C _Peq] is then given by Eq. (2.16), with the association rate constant k _a replaced by the activation-controlled limit k _A. In Laplace space, the relaxation function is

(4.59)

The second limiting situation is when binding and unbinding are under the influence of diffusion but occur infrequently. Then a newly released ligand molecule will have sufficient time, before the next binding event, to equilibrate with the other ligand molecules and reach the steady-state distribution of irreversible binding. Again, the ordinary rate equation applies for the relaxation of the protein concentration at times longer than that for reaching the steady-state distribution of irreversible binding, but now the association rate constant is the k _a calculated in the preceding subsections under the assumption of irreversible binding.

When binding and unbinding are under the influence of diffusion and not infrequent, a newly released ligand molecule and the other ligand molecules do not have sufficient time to relax to the steady-state distribution of irreversible binding; the concentration of ligand molecules around the protein will thus be somewhat higher than expected from that distribution, and hence the effective binding rate constant will be higher than the k _a for irreversible binding.

Gopich & Szabo (Reference Gopich and Szabo2002) devised a self-consistent relaxation time approximation for the kinetics of reversible diffusion-influenced binding. In Laplace space, their relaxation function is given by

(4.60a)

where

(4.60b)

has a structure similar to a gated rate coefficient for irreversible binding [see Eq. (4.49)]. The ‘gating’ rate k ₀ is determined self-consistently by requiring that the area under the relaxation function, i.e., Ŷ(0), is the inverse of k ₀. When diffusion is infinitely fast, s _a(s) and s _G(s)→k _A and Eq. (4.60a) reduces to Eq. (4.59). The connection to a gated rate coefficient for irreversible binding comes from the fact that the unbound protein allows the surrounding ligand molecules to bind but the bound protein is reflecting to them. Notice, however, that Eq. (4.60b) differs from Eq. (4.49) by the appearance of ρ_beq on the left-hand side. As a result, here s _G(s)>s _a(s), an outcome already anticipated in the preceding paragraph. In contrast, gating transition from an active form to an inactive from always reduces the rate coefficient for irreversible binding.

5.1 Effects on diffusion constants

It seems that not much attention has been paid by theorists to how internal dynamics of macromolecules is affected by crowding, although the problem has been studied by molecular dynamics simulations (Cheung et al. Reference Cheung, Klimov and Thirumalai2005; Minh et al. Reference Minh, Chang, Trylska, Tozzini and Mccammon2006; Mittal & Best, Reference Mittal and Best2010; Qin et al. Reference Qin, Minh, Mccammon and Zhou2010). On the other hand, the problem of how the translational diffusion constant of a macromolecule is affected by crowding is a classical problem. The theory of Tokuyama & Oppenheim (Reference Tokuyama and Oppenheim1994, Reference Tokuyama and Oppenheim1995), for a tracer hard sphere crowded by other identical hard spheres, is one of the most rigorous. Each particle undergoes Brownian motion, with diffusion constant D, were it not for hydrodynamic interactions and collisions with other particles. The effective diffusion constant of the tracer particle, as a function of the crowder volume fraction Φ, is given by

(5.1a)

where

(5.1b)

(5.1c)

with b ₁=(9Φ/8)^1/2, b ₂=11Φ/6, and Φ₀≈0·5718. The numerator of the right-hand side of Eq. (5.1a) accounts for collisions, the L _l(Φ) term accounts for local hydrodynamic interactions (arising from a static distribution of neighboring particles), and the L _n(Φ) term accounts for non-local hydrodynamic interactions (arising from a relaxed distribution of neighboring particles). For a mixture of different-sized ‘soft’ crowders, Tokuyama (Reference Tokuyama2009a, Reference Tokuyamab) suggests fitting experimental data according to the following formula:

(5.2)

with η and Φ₀ as adjustable parameters.

At Φ=0·3, which approximates the level of intracellular crowding (Zimmerman & Trach, Reference Zimmerman and Trach1991), Eq. (5.1a) predicts a 5-fold reduction in diffusion constant. This is generally consistent with the extent of crowding-induced reduction in translational diffusion constant measured on proteins in cytoplasm and concentrated solutions (Kuttner et al. Reference Kuttner, Kozer, Segal, Schreiber and Haran2005; Li et al. Reference Li, Wang and Pielak2009; Swaminathan et al. Reference Swaminathan, Hoang and Verkman1997). These experiments also show that rotational diffusion constants are reduced to a lesser extent, as to be expected. Similarly, one expects the effect of crowding on the effective friction coefficient or diffusion constant for internal dynamics to be relatively modest.

5.2 Effects on free energies of transition states and transient complexes

Crowding also affects the energy functions of macromolecules. A universal and significant aspect of the crowding effect, known as excluded-volume interactions, arises from the fact that molecules cannot occupy the same region in space. Excluded-volume interactions increase the free energy of a macromolecule, but the increase depends on its conformation. Compared to a more open conformation, a compact conformation is more easily accommodated into a distribution of crowders, and hence the crowding-induced increase in free energy, ΔG _c, is less. Effectively, crowding changes the relative stability between these two conformations in favor of the latter. According to scaled particle theory (Lebowitz et al. Reference Lebowitz, Helfand and Praestgaard1965), ΔG _c for a spherical test particle (with radius R _p) crowded by other spherical crowders (with radius R _c and occupying a volume fraction Φ) is

(5.3)

where C is the number density of the crowders, Σ_c and Σ_p are the surface areas of the crowders and the test particle, respectively, and V _p is the volume of the test particle. As to be expected, ΔG _c increases with increasing R _p. A recent development extends Eq. (5.3) to atomistic proteins, with the geometric parameters R _p, Σ_p, and V _p calculated from protein conformations (Qin and Zhou, Reference Qin and Zhou2010).

When a protein undergoes a folding transition, the transition state is more compact than the unfolded state. Excluded-volume interactions with crowders are thus expected to stabilize the folding transition state. If it is assumed that the effective diffusion constant along the folding reaction coordinate is not affected (see preceding subsection), then crowding would increase the folding rate.

Similarly, crowding is expected to stabilize the transient complex of a protein–protein pair relative to the unbound state. Crowding thus produces two opposing contributions to the diffusion-controlled association rate constant: a decrease due to the reduction in translational (and to a less extent rotational) diffusion constant and an increase due to the stabilization of the transient complex. The two opposing contributions lead to a rather modest net effect by crowding. Assuming that the crowding-induced change in free energy, ΔΔG _c, with the unbound state as a reference, is long ranged, one may adapt Eq. (4.46) to obtain the following result for the association rate constant under crowding (Zhou, Reference Zhou2004):

(5.4)

where k _a is the rate constant in a dilute solution and <ΔΔG _c>* is the average of ΔΔG _c in the transient complex. [An additional factor, arising from the disparate effects of crowding on translational and rotational diffusion constants, may also be introduced; see Eq. (4.36d).] Calculations suggest that ΔΔG _c, like electrostatic attraction between two proteins, is indeed long ranged (Qin, Lu, and Zhou, to be published).

6.1 Unfolding/unzipping under force

The unfolding of proteins, unzipping of nucleic acids, and dissociation of protein–ligand complexes have been studied in many single-molecule ‘pulling’ experiments (Cao et al. Reference Cao, Kuske and Li2008; Carrion-Vazquez et al. Reference Carrion-Vazquez, Oberhauser, Fowler, Marszalek, Broedel, Clarke and Fernandez1999; Chen et al. Reference Chen, Wen and Tinoco2007; Greenleaf et al. Reference Greenleaf, Frieda, Foster, Woodside and Block2008; Kellermayer et al. Reference Kellermayer, Smith, Granzier and Bustamante1997; Liphardt et al. Reference Liphardt, Onoa, Smith, Tinoco and Bustamante2001; Marshall et al. Reference Marshall, Sarangapani, Lou, Mcever and Zhu2005; Mathé et al. Reference Mathé, Visram, Viasnoff, Rabin and Meller2004; Merkel et al. Reference Merkel, Nassoy, Leung, Ritchie and Evans1999; Rief et al. Reference Rief, Gautel, Oesterhelt, Fernandez and Gaub1997, Reference Rief, Clausen-Schaumann and Gaub1999; Schlierf & Rief, Reference Schlierf and Rief2006; Yang et al. Reference Yang, Cecconi, Baase, Vetter, Breyer, Haack, Matthews, Dahlquist and Bustamante2000). As the pulling force is increased monotonically with time, at some point the system undergoes an unfolding/unzipping (or dissociation) transition. The pulling force at that point is referred to as the rupture force. The transition can be modeled as a rate process, with the rate coefficient k _u affected by the time-dependent pulling force F(t). The survival probability of the intact system is then governed by (Evans & Ritchie, Reference Evans and Ritchie1997)

(6.1a)

which has the solution

(6.1b)

For a stiff system pulled at a constant speed v by a spring with spring constant κ_s, one has F=κ_svt and =κ_sv. The probability density of the rupture force, ρ(F), is related to the survival probability via

(6.2a)

Hence,

(6.2b)

If k _u(F) is calculated from a theory (see below), Eq. (6.2b) can be used to predict the distribution of the rupture force.

Experimentalists face the reverse problem: they measure the distribution of the rupture force and want to extract k _u(F) from their data.Footnote ⁵ To that end, one rearranges Eq. (6.1a) to express k _u(F) in terms of S(t) and its time derivative and then uses Eq. (6.2a) to express the latter two quantities in terms of ρ(F). The result is

(6.3)

As observed by Dudko et al. (Reference Dudko, Hummer and Szabo2006, Reference Dudko, Hummer and Szabo2008), data collected at different pulling speeds and analyzed according to Eq. (6.3) should all collapse to the same curve.

Now we consider the calculation of k _u(F). The unfolding/unzipping transition can be modeled as an escape from a potential well along the pulling direction. The pulling force modifies the potential function from U ₀(x) to U ₀(x)−Fx, where x is the coordinate along the pulling direction. Assuming that the motion along x is diffusive, Kramers' formula [Eq. (3.44)] predicts the force-dependent rate coefficient as (Evans & Ritchie, Reference Evans and Ritchie1997; Shapiro & Qian, Reference Shapiro and Qian1997)

(6.4)

where k _u0 is the rate coefficient in the absence of force, and well and barrier mean that the integrations are restricted to these regions. Dudko et al. (Reference Dudko, Hummer and Szabo2006) evaluated these integrals approximately for a quadratic-cusp form and a linear-cubic form of U ₀(x). The results in both cases can be written as

(6.5)

where ν=1/2 and 2/3 for the two forms of potentials, respectively; ΔU‡ is the activation energy of U ₀(x); and Δx‡ is the displacement between the bottom of the well and the top of the barrier. When ν=1, Eq. (6.5) reduces to Bell's formula (Bell, Reference Bell1978).

Dudko et al. (Reference Dudko, Hummer and Szabo2008) analyzed data from two single-molecule experiments. One studied the unfolding of a protein by atomic force microscopy (Schlierf & Rief, Reference Schlierf and Rief2006); the other studied the voltage-induced unzipping of a DNA pair in a nanopore (Mathé et al. Reference Mathé, Visram, Viasnoff, Rabin and Meller2004). The data from each experiment, collected at different pulling or voltage-ramp speeds, when processed according to Eq. (6.3), collapsed to a single curve, yielding k _u(F). They then fitted these k _u(F) to Eq. (6.5) to generate values for the microscopic parameters k _u0, ΔU‡, and Δx‡.

6.2 Kinetics of fluctuating enzymes

Single-molecule experiments have revealed that a single enzyme can produce a wide range of turnover rates, presumably due to conformational fluctuations (Antikainen et al. Reference Antikainen, Smiley, Benkovic and Hammes2005; English et al. Reference English, Min, Van Oijen, Lee, Luo, Sun, Cherayil, Kou and Xie2006; Flomenbom et al. Reference Flomenbom, Velonia, Loos, Masuo, Cotlet, Engelborghs, Hofkens, Rowan, Nolte, Van Der Auweraer, De Schryver and Klafter2005; Lu et al. Reference Lu, Xun and Xie1998; Shi et al. Reference Shi, Dertouzos, Gafni, Steel and Palfey2006; van Oijen et al. Reference Van't Hoff2003; Zhang et al. Reference Zhang, Rajagopalan, Selzer, Benkovic and Hammes2004). In particular, English et al. (Reference English, Min, Van Oijen, Lee, Luo, Sun, Cherayil, Kou and Xie2006) followed events of turning over a fluorogenic substrate by β-galactosidase. The turnover rates at each substrate concentration fluctuate widely, but the dependence of the average rate, as measured by the inverse of the average waiting time, 〈τ〉, between turnover events, on substrate concentration conforms to the Michaelis–Menten relation. This apparent Michaelis–Menten behavior initially was explained by invoking the assumption that conformational exchanges in the enzyme–substrate complex are extremely slow (i.e., in the quasi-static limit) (Kou et al. Reference Kou, Cherayil, Min, English and Xie2005). A later reexamination suggested that the Michaelis–Menten behavior was the result of a very small k _cat (i.e., satisfying the quasi-equilibrium condition) (Min et al. Reference Min, Gopich, English, Kou, Xie and Szabo2006; see also Qian, Reference Qian2008).

The sequence of waiting times contains rich information. Besides the average, the distribution function φ(τ) and the correlation function c_i=〈(τ_i−〈τ〉)(τ₀−〈τ〉)〉, where τ_i is the ith waiting time, are useful. Min & Xie (Reference Min and Xie2006) showed that the correlation function measured on β-galactosidase (English et al. Reference English, Min, Van Oijen, Lee, Luo, Sun, Cherayil, Kou and Xie2006) is mimicked by the correlation between events of crossing a tiny (0·1 k _BT) energy barrier in a double-well potential, when the dynamics is governed by the generalized Langevin equation with a power-law memory kernel. That memory kernel itself was motivated by the observed dynamic behavior of the distance between the electron donor and acceptor in a protein (Min et al. Reference Min, Luo, Cherayil, Kou and Xie2005). The fluctuations of the distance were modeled well as a particle moving in a harmonic well according to the generalized Langevin equation with a power-law memory kernel.

6.3 Protein folding and association under crowding

Ai et al. (Reference Ai, Zhou, Bai and Choy2006) measured the folding and unfolding rate constants of a four-helix bundle protein under crowding. In the presence of 85 g l⁻¹ of a crowding agent, PEG 20 K, the folding rate constant is increased by 2-fold but the unfolding rate constant is essentially unchanged. Qualitatively, these observations are consistent with crowding effects dominated by excluded-volume interactions. To carry out a quantitative test, the transition-state ensemble of this protein was generated from molecular dynamics simulations (Tjong & Zhou, Reference Tjong and Zhou2010). It consists of highly compact conformations with residual secondary structures. Calculations of crowding effects on the transition-state ensemble yield changes in folding and unfolding rate constants that are in agreement with experimental results.

Yuan et al. (Reference Yuan, Chyan, Zhou, Chung, Peng, Ping and Yang2008) studied the effects of crowding on the pulling force required to unfold ubiquitin. Compared to the result in a dilute solution (unfolding force at 210 pN), the unfolding force increases by 24 pN when the concentration of a crowding agent, dextran 40 K, is increased to 300 g l⁻¹. The increase in pulling force was interpreted as the reflection of a decrease in the unfolding rate constant. The latter in turn was attributed to the destabilization of the transition state, relative to the more compact folded state, by dextran crowding. By representing the transition state as a sphere somewhat enlarged from the folded state and using Eq. (5.3) to calculate the change in activation energy by crowding, the experimental results could be explained.

Kuttner et al. (Reference Kuttner, Kozer, Segal, Schreiber and Haran2005) studied how PEG 8 K affects the association rate constant of two proteins. To dissect the contributions of crowding, they also measured the effects of the crowding agent on translational and rotational diffusion constants of proteins. As already noted in subsection 5.1, crowding affects the two types of diffusion constants to different extents. After correcting for the reductions in diffusion constants, Kuttner et al. found that the association constant is higher than expected, implicating a crowding-induced stabilization of the transient complex; the stabilization is ~0·8 k _BT, in line with calculations (Qin, Lu, and Zhou, to be published).