format: text

text:
# Microscopic Theory of Protein Folding Rates.II: Local Reaction Coordinates and Chain Dynamics

## I Introduction

A complete description of protein folding kinetics must account for the microscopic dynamics of the polypeptide chain. The relaxation of a free polymer is well described by Langevin dynamics with an harmonic potential between adjacent monomers and Markovian random forces . This Rouse-Zimm formalism can be extended through frictional forces with memory to effectively describe chain motions that include many other strictly microscopic complexities such as chain stiffness . To describe the chain motions involved in protein folding, one must also attend to the residue specific interactions or free energy landscape which guides the protein to its native state and governs the probability of different ensembles of structures described in terms of order parameters.

For many proteins, the only thermodynamically distinct states with significant detectable populations are the unfolded ensemble and the native folded conformations . This observation suggests the existence of a free energy barrier separating the stable and metastable phases. Like nucleation in phase transitions, the folding kinetics for these proteins is largely governed by the probability of rare configurations corresponding to the top of the barrier, the so-called transition state ensemble. In theoretical studies, the transition state ensembles correspond with saddle-points in a multi-dimensional free energy surface as a function of local order parameters. Near the transition state, some order parameter fluctuations grow to become the folded phase. This is an unstable mode of collective motion. Motion in other directions tend to stabilize the transition state. The structure of the unstable growth mode can be thought of as a locally defined reaction coordinate. If there are multiple saddle points which are connected, the local reaction coordinate for each of these will usually differ. In such cases, especially, the local reaction coordinate may not be the best way of describing the global transformation of the ensemble involved in folding.

Lately, there has been a great effort to elucidate the nature of the transition state ensemble and ideal reaction coordinate for specific proteins through measurement , theoretical calculations based on physically intuitive choices , and simulations . While many aspects of the folding kinetics seem to be understood, still some basic issues remain controversial. For example, while the structural distribution of the transition state ensemble is well predicted , the absolute magnitude of the barrier height is not well determined. Comparison with experiment is clouded by the uncertainty of the prefactor of the folding rate.

For simple thermally activated transitions, one expects the rate to follow the activated law

$$k=k_0e^{\beta \mathrm{\Delta }F^{}},$$
(1)
where $`\beta =1/k_BT`$ is the inverse temperature and $`\mathrm{\Delta }F^{}`$ is the barrier height which contains both entropy and energy. For folding, $`\mathrm{\Delta }F^{}`$ is expected to be strongly temperature dependent. The exponential factor is the relative probability of finding the system at the transition state compared to the metastable minima. The prefactor, $`k_0`$, is determined by the microscopic dynamics involved in crossing the barrier. For proteins with very rugged landscapes, the dynamics can be described as hopping between trapped states . In this case the barriers to escape traps enter strongly. For proteins with very little frustration that have smoother landscapes, there will be little trapping, but the chain motion will be diffusive and depend on chain connectivity. $`k_0`$ depends on the free energy surface and the frictional forces such as viscosity that determine the stochastic motion of the protein as it folds. One estimate proposed for the prefactor relating $`k_0`$ to the inverse time it takes to form a typical loop given by simple polymer theory, $`k_01\mu \mathrm{s}^1`$ . While quite reasonable, this estimate does not involve the shape of the free energy surface near the transition state. Thus, for example, it is not completely clear how big a loop should be considered.

In Ref. jjp:st:pgw:98  and the preceding paper (henceforth referred as \[I\]), we presented a variational theory for the protein folding free energy surface of a specific protein to a known structure. In this approach, the protein is modeled as a collapsed, stiff chain with interactions between residues in contact in the native structure. The reference Hamiltonian is that of a polymer chain inhomogeneously constrained to the native structure. The free energy profile for folding obtained from model can be described as an activation barrier between the unfolded and folded ensembles with fine structure from local minima and saddle-points in the free energy surface giving rise to longitudinal ruggedness in the free energy profile for a perfectly funneled landscape. In this approach, conformational ensembles along the folding route are characterized by structural order parameters evaluated as averages over the distribution from the reference Hamiltonian. The focus of the present paper is how the microscopic chain dynamics are involved in barrier crossing within this model.

Understanding barrier crossing dynamics in protein folding requires not only a good characterization of the free energy surface and structural ensembles, but also a good description of the microscopic motions of the polymer backbone. In the Rouse-Zimm model of polymer dynamics, the motion of the polymer chain is governed by the chain connectivity, represented as harmonic entropic springs connecting adjacent monomers, and is damped by the interactions of the monomers with the solvent . A more detailed description of polypeptide dynamics acknowledges there are activated transitions between gauche and trans bond angles determined by the dihedral potential . One way to incorporate the short range structural information arising from these barriers to internal rotation is through a static chain stiffness that defines renormalized harmonic correlations . However, a complete description must also incorporate the activated nature of these rotational isomerization rates in the dynamical correlations.

Quite early in the study of polymer dynamics, anomalies found in intrinsic viscosity measurements indicated that there was a source of friction independent of the viscosity of the solvent . Kuhn and Kuhn first suggested this “internal viscosity” was a consequence of microscopic barrier hopping in the local chain segments . This phenomenological friction found for high polymers may have a variety of microscopic origins in addition to internal friction due to barrier hopping of many elementary chain links. A helpful (though perhaps dated) evaluation of different theories can be found found in Ref. cwm:mcw:85 . Fixman has shown that simulations for polymers with rigidly fixed bond angles and bond lengths can be mimicked by internal viscosity . Still other possible sources of the measured internal viscosity include the effect of interactions with distant monomers in the chain , or other inherent non-linearities neglected in the Gaussian model assumption . To our knowledge, the microscopic interpretation of internal viscosity for artificial high polymers has not been definitively settled. It is clear, however, that polypeptides will at least exhibit internal viscosity from the difficulty of making backbone angle transitions. These transitions are roughly 100 times slower than if no barriers were present in the Ramachandran potential.

In this paper, we complete the logical development of folding rate theory for minimally frustrated proteins by analyzing the microscopic dynamics of chain motions involved in barrier crossing of the protein folding model presented in Ref. jjp:st:pgw:98  and \[I\]. We use multi-dimensional, non-Markovian Langevin dynamics to describe the barrier crossing, and identify the unstable modes as the reaction coordinates. For the folding rate calculation, we use a generalization of Kramer’s theory  that treats non-Markovian dynamics in many dimensions . At the same time, the interpretation of a local reaction coordinate follows Langer’s theory of the nucleation rates . The details of chain dynamics enter the formalism through a generalized friction matrix with memory. The friction that damps the motion of the chain includes the solvent viscosity and an internal friction that depends on the length scale (in sequence) of the motion but which we take as is independent of the solvent viscosity. The internal friction accounts for the activated nature of dihedral angle flips in polypeptides. We restrict our attention to over-damped dynamics, though it is easy to generalize the calculations to include inertial terms as well.

## II Reaction Coordinates and Folding Rates

To calculate the barrier crossing rate, one must first identify the order parameters that describe the phases of the system, and then study their motions. In density functional theory, one often considers a density field that is slowly varying in space, and expands the density functional about a uniform density . This expansion yields a free energy consisting of a bulk free energy and can be expanded in gradient terms that ultimately represent the surface tension between two phases. The density then evolves either hydrodynamically or diffusively under the influence of the free energy driving force. The quasi-equilibrium part of this approach has been applied to nucleation of the folded state in proteins starting with a mean field free energy functional and defining a global native density field as a measure of similarity to the native state

$$\rho ^\mathrm{N}(𝐫,\{𝐫_i\})\underset{i}{}\delta (𝐫𝐫_i)\delta (𝐫_𝐢𝐫_i^\mathrm{N}),$$
(2)
where $`\{𝐫_i\}`$ specify the configuration of the protein chain, and $`\{𝐫_i^\mathrm{N}\}`$ are the corresponding coordinates of the native structure. In the present work, we wish to have a description of the folding that is local in sequence and hence consider each term in the sum separately. This local native density defined as

$$\rho ^\mathrm{N}(𝐫,𝐫_i)=\delta (𝐫𝐫_i)\mathrm{exp}\left[\frac{3}{2a^2}\alpha ^\mathrm{N}(𝐫_i𝐫_i^\mathrm{N})^2\right],$$
(3)
where we have relaxed the measure of similarity of the native state by replacing $`\delta (𝐫_𝐢𝐫^\mathrm{N})`$ by a Gaussian measure.

To connect this density to our folding route calculations, it is more convenient to consider scalar order parameters than to expand the free energy in terms of field variables. Therefore, we integrate $`\rho (𝐫,𝐫_i)`$ over $`𝐫`$, to give the scalar measure native similarity

$$\rho (𝐫_i)=\mathrm{exp}[\frac{3}{2a^2}\alpha ^\mathrm{N}(𝐫_i𝐫_i^\mathrm{N})^2]].$$
(4)
\[To simplify notation, we have suppressed the superscript $`N`$ since the native density can still be distinguished by the 1-body monomer density, $`\rho ^1(𝐫)`$, given in Eq.(I-25)\]. Following the interpretation of the order parameters given in \[I\], we calculate the average native density of site $`i`$, $`\rho _i=\rho (𝐫_i)_0`$, for a given set of variational constraints $`\{C_i\}`$

$`\rho _i[\{C\}]={\displaystyle 𝑑𝐫\rho _i^1(𝐫)\rho (𝐫)}`$ (6)
$`=(1+\alpha ^\mathrm{N}G_{ii})^{3/2}\mathrm{exp}\left[{\displaystyle \frac{3}{2a^2}}{\displaystyle \frac{\alpha ^\mathrm{N}(𝐬_i𝐫_i^\mathrm{N})^2}{1+\alpha ^\mathrm{N}G_{ii}}}\right];`$
that is, $`\rho _i`$ is a function of $`\{C_i\}`$ through the correlations, $`G_{ij}`$, and average position, $`𝐬_i`$, given by Eq.(I-21) and Eq.(I-22), respectively. The free energy can then be parameterized by $`\{\rho _i\}`$ with $`\overline{F}[\{\rho \}]=F[\{C\}]`$ where $`\{\rho \}`$ is evaluated at $`\{C\}`$. In particular, we denote a saddle-point of $`F[\{C\}]`$ by $`\{C^{}\}`$ and the corresponding native density by $`\rho _i^{}=\rho _i[\{C^{}\}]`$.

We assume Langevin dynamics with memory for the native density, $`\{\rho _i(t)\}`$,

$$_t\rho _i(t)=\underset{j}{}_0^t𝑑t^{}\mu _{ij}(tt^{})\frac{\beta \overline{F}[\{\rho (t^{})\}]}{\rho _j}+\xi _i(t),$$
(7)
where $`\beta =1/k_BT`$ is the inverse temperature, and $`\mu (t)`$ is a generalized mobility matrix related to the random noise $`\xi _i(t)`$ through the correlations $`\xi _i(t)\xi _j(t^{})=\mu _{ij}(tt^{})`$. The form of $`\mu (t)`$ is discussed in the following section, here we assume that $`\mu (t)`$ is a known function of time.

To study the dynamics near the saddle-point, we expand $`\overline{F}[\{\rho \}]`$ about $`\rho _i^{}`$ to second order

$$\beta \overline{F}[\{\rho \}]\mathrm{const}+\frac{1}{2}\underset{ij}{}\overline{\mathrm{\Gamma }}_{ij}^{}\delta \rho _i\delta \rho _j,$$
(8)
where $`\delta \rho _i=\rho _i\rho _i^{}`$, and the Hessian matrix

$$\overline{\mathrm{\Gamma }}_{ij}^{}=\frac{^2\beta \overline{F}[\{\rho ^{}\}]}{\rho _i\rho _j}$$
(9)
has one negative eigenvalue since it is evaluated at a saddle-point. In terms of derivatives with respect to the variational parameters, $`\overline{\mathrm{\Gamma }}^{}`$ is the solution of the matrix equation

$$\underset{kl}{}\overline{\mathrm{\Gamma }}_{kl}^{}\frac{\rho _k^{}}{C_i}\frac{\rho _l^{}}{C_j}=\frac{^2\beta F[\{C^{}\}]}{C_iC_j}.$$
(10)
where the derivatives with respect to $`C_i`$ can calculated as described in \[I\].

With this approximation, Eq.(7) can be written as

$$_t(\delta \rho _i(t))=\underset{j}{}_0^t𝑑t^{}\left[\mu (tt^{})\overline{\mathrm{\Gamma }}^{}\right]_{ij}\delta \rho _j(t^{})+\xi _i(t).$$
(11)
Accordingly, the native density correlation functions $`𝒞(t)`$ with components

$$𝒞_{ij}(t)=\delta \rho _i(t)\delta \rho _j(0)_0,$$
(12)
satisfy

$$_t𝒞(t)=_0^t𝑑t^{}\mu (tt^{})\overline{\mathrm{\Gamma }}^{}𝒞(t^{}).$$
(13)
In Eq.(11), we have subtracted the stationary value $`\rho _i(t\mathrm{})=\rho _i^{}`$, so that $`\delta \rho _i(t\mathrm{})=0`$ and $`𝒞(t\mathrm{})=0`$. This long time solution is unstable, however, because it is a saddle-point of the free energy.

The unstable mode can be determined by the average equation of motion . Denoting the Laplace transform of an arbitrary function of time by $`\widehat{g}(\omega )=_0^{\mathrm{}}𝑑te^{\omega t}g(t)`$, the solution of the Eq.(11) averaged over the random noise is

$$\overline{\delta \widehat{\rho }_i(\omega )}=\left[\omega \mathrm{𝟏}+\widehat{\mu }(\omega )\overline{\mathrm{\Gamma }}^{}\right]_{ij}^1\delta \rho _j(0).$$
(14)
By inverting the Laplace transform, $`\overline{\rho _i(t)}`$ can be expressed as a sum of exponentials with time constants determined by the poles of $`\left[\omega \mathrm{𝟏}+\widehat{\mu }(\omega )\overline{\mathrm{\Gamma }}^{}\right]^1`$: the eigenfunctions $`\widehat{\mu }(\kappa )\overline{\mathrm{\Gamma }}^{}𝐮=\kappa 𝐮`$ have time dependence $`𝐮(t)e^{\kappa t}`$. Assuming that $`\widehat{\mu }(\omega )`$ is positive definite, there is one mode for which $`\kappa `$ is negative, because $`\overline{\mathrm{\Gamma }}`$ is the saddle point Hessian of $`F[\rho ]`$:

$$𝐮^{}(t)=𝐮^{}e^{\left|\kappa ^{}\right|t},$$
(15)
where

$$\widehat{\mu }(\kappa ^{})\overline{\mathrm{\Gamma }}^{}𝐮_j^{}=\kappa ^{}𝐮_j^{}(\kappa ^{}<0).$$
(16)
Since motion along the unstable mode $`𝐮^{}`$ grows exponentially away from the saddle point, we identify the components of $`𝐮^{}`$ as the local reaction coordinate to surmount the barrier. Both $`𝐮^{}`$ and $`\kappa ^{}`$ are essentially the unstable mode and curvature of $`\overline{F}[\rho ]`$, but renormalized by the dynamics of the barrier crossing incorporated in the mobility matrix. In Eq.(16), $`\kappa ^{}`$ is the multi-dimensional analogue of the Grote-Hynes frequency (in the over-damped limit) .

The rate for barrier crossing corresponding to Eq.(11) is

$$k=\frac{|\kappa ^{}|}{2\pi }\left|\frac{det\overline{\mathrm{\Gamma }}_{MS}}{det\overline{\mathrm{\Gamma }}^{}}\right|^{1/2}e^{\beta \mathrm{\Delta }F^{}},$$
(17)
where $`\mathrm{\Delta }F^{}`$ is the barrier height and $`\overline{\mathrm{\Gamma }}_{MS}`$ and $`\overline{\mathrm{\Gamma }}^{}`$ are the curvature matrices of $`\overline{F[\rho ]}`$ evaluated at the metastable minimum and the saddle-point, respectively. Eq.(17) a generalizes both the rate calculations presented by Langer  and Grote-Hynes ; it simultaneously accommodates both multi-dimensional diffusion (as in the Langer formula) as well as time dependent friction (as in the Grote-Hynes formula).

The ratio of determinants accounts for the entropic differences between the transition state and the metastable phase due to fluctuations of the order parameter. These fluctuations are like the capillary waves which renormalize the surface tension in ordinary nucleation. These contributions would already be included in an exact free energy so that the determinants should be absorbed into the exponential factor to most simply keep a consistent level of thermodynamic theory . Consequently, the expression for the rate becomes

$$k=\frac{|\kappa ^{}|}{2\pi }e^{\beta \mathrm{\Delta }F^{}},$$
(18)
where $`\kappa ^{}`$ is given by Eq.(16). As pointed out in Ref. dfc:pgw:83 , if the over-counting of entropy were not corrected, the ratio of the forward and backwards rate would not equal the equilibrium constant as determined by the starting free energy functional, $`k_{12}/k_{21}e^{\beta (F_2F_1)}`$.

## III Linear Response Approximation of the Generalized Mobility

The time evolution of the native densities $`\{\delta \rho _i(t)\}`$ depends on microscopic dynamics of the monomer positions $`\{𝐫_i(t)\}`$. The chain dynamics, in turn, are incorporated into the effective mobility matrix $`\mu (t)`$ in the relaxation equation for $`\delta \rho (t)`$ (Eq.(11)). To determine $`\mu (t)`$ consistent with the definition of the native densities this theory, we consider the chain dynamics near the saddle point to be governed by the $`H_0`$ with constraints $`\{C^{}\}`$. The resulting mobility matrix for native densities in Eq.(11) is assumed to be same as that which arises from the Mori-Zwanzig projected dynamics . We begin by outlining this formalism as it was applied in Zwanzig’s classic paper on generalized Rouse dynamics for a polymer described by an arbitrary potential . A closed form of the mobility matrix is then obtained by approximating the local dynamics with those of the constrained polymer Hamiltonian, $`H_0`$.

Consider a polymer described by the potential $`U(𝐑)`$ with monomer positions $`𝐑=(𝐫_1,\mathrm{},𝐫_n)`$. The probability density of monomer positions, $`\mathrm{\Psi }(𝐑,t)`$ is assumed to evolve according to the Smoluchowski equation

$$_t\mathrm{\Psi }(𝐑,t)=𝒟\mathrm{\Psi }(𝐑,t)$$
(19)
$$𝒟=\underset{ij}{}_iD_{ij}e^{\beta U(𝐑)}_je^{\beta U(𝐑)}$$
(20)
where $`_i/𝐫_i`$, and $`D`$ is the diffusion matrix. In this section, we denote the equilibrium averages by $`\mathrm{}=𝑑𝐑\mathrm{}\mathrm{\Psi }_{eq}(𝐑)`$ where

$$\mathrm{\Psi }_{\mathrm{eq}}(𝐑)=e^{\beta U(𝐑)}\left[𝑑𝐑e^{\beta U(𝐑)}\right]^1.$$
(21)
Now consider a set of dynamical variables that are functions of the monomer positions, $`\{A_i(𝐑)\}`$. These variables will be chosen to be the order parameters $`\{\rho \}`$ to describe the barrier crossing dynamics discussed in the previous section, but for now $`\{A\}`$ is arbitrary. The equilibrium time dependent correlation functions of

$$\delta A_i=A_iA_i$$
(22)
are denoted by

$$𝒞_{ij}^A(t)\delta A_i(t)\delta A_j(0).$$
(23)
The time dependent correlations can be formally expressed as

$`𝒞_{ij}^A(t)`$ $`=`$ $`{\displaystyle 𝑑𝐑\delta A_ie^{𝒟t}(\delta A_j\mathrm{\Psi }_{\mathrm{eq}})}`$ (24)
$`=`$ $`{\displaystyle 𝑑𝐑\mathrm{\Psi }_{\mathrm{eq}}\delta A_ie^t\delta A_j}.`$ (25)
In the second line of Eq(24), we have introduced the adjoint operator $``$ defined by its action on an arbitrary function $`B`$

$$𝒟\mathrm{\Psi }_{\mathrm{eq}}B=\mathrm{\Psi }_{\mathrm{eq}}B,$$
(26)
or explicitly,

$$=\underset{ij}{}e^{\beta U(𝐑)}_iD_{ij}e^{\beta U(𝐑)}_j.$$
(27)
The meaning of the operator notation in Eq.(24) is defined through the propagator (or Green’s function) $`P(𝐗t;𝐘)`$:

$$𝒞_{ij}^A(t)=𝑑𝐗𝑑𝐘\delta A_i(𝐗)\delta A_j(𝐘)P(𝐗t;𝐘)\mathrm{\Psi }_{\mathrm{eq}}(𝐘),$$
(28)
where (formally)

$$P(𝐗t;𝐘)=e^{𝒟t}\delta (𝐗𝐘),$$
(29)
i.e., $`P(𝐗t;𝐘)`$ satisfies the Smoluchowski equation, $`_tP=𝒟P`$, with initial conditions $`P(𝐗t=0;𝐘)=\delta (𝐗𝐘)`$.

An alternative equation can be derived for the correlations through the projected dynamics of Eq.(19). Following Ref. zwanzig:74 , we define the projection operator through its action on a arbitrary function $`B`$

$$𝒫B=\underset{ij}{}\delta A_i\mathrm{\Gamma }_{ij}^\mathrm{A}\delta A_jB$$
(30)
where $`(\mathrm{\Gamma }^\mathrm{A})^1`$ is the matrix of static correlations

$$(\mathrm{\Gamma }^\mathrm{A})^1=\delta A_i\delta A_j.$$
(31)
Using standard projection operator techniques , it can be shown that the correlations obey the generalized Langevin equation

$$_t𝒞^\mathrm{A}(t)=\mathrm{\Omega }^\mathrm{A}\mathrm{\Gamma }^\mathrm{A}𝒞^\mathrm{A}+_0^t𝑑t^{}K^\mathrm{A}(tt^{})\mathrm{\Gamma }^\mathrm{A}𝒞^\mathrm{A}(t^{}),$$
(32)
where

$$\mathrm{\Omega }_{ij}^\mathrm{A}=A_iA_j=\underset{kl}{}_kA_iD_{kl}_lA_j$$
(33)
and $`K(t)`$ is the memory kernel

$$K_{ij}^\mathrm{A}(t)=A_i(1𝒫)\mathrm{exp}[(1𝒫)t](1𝒫)A_j.$$
(34)
Eq.(32) can be written in the form of Eq.(12)

$$_t𝒞^\mathrm{A}(t)=_0^t𝑑t^{}\mu ^\mathrm{A}(tt^{})\mathrm{\Gamma }^\mathrm{A}𝒞^\mathrm{A}(t^{}),$$
(35)
where

$$\mu ^\mathrm{A}(t)=2\mathrm{\Omega }^\mathrm{A}\delta (t)K^\mathrm{A}(t)$$
(36)
is the effective time-dependent mobility matrix. As suggested by this notation, the static correlations $`(\mathrm{\Gamma }^\mathrm{A})^1`$ can be thought of as arising from an effective harmonic free energy $`\beta F=\mathrm{const}+(1/2)_{ij}\delta A_i\mathrm{\Gamma }_{ij}^\mathrm{A}\delta A_j`$ .

Eq.(35) is an exact result, equivalent to Eq.(28). The linear response formulation avoids explicit use of the Green’s function which is generally unknown for an arbitrary potential. On the other hand, it is very difficult to calculate the memory kernel explicitly from Eq.(34) since the projection operators are very awkward to manipulate.

If we choose $`\{A\}`$ to be the monomer positions, then Eq.(32) becomes an equation for the monomer correlation. If one neglects the memory term, the monomer correlations are said to obey optimized Rouse dynamics . This is equivalent to having an effective pre-averaged diffusion matrix, $`\mathrm{\Omega }=D`$, and a Gaussian approximation to the chain through the static correlations. \[Bixon and Zwanzig first presented the stiff chain connectivity matrix \[Eq.(I-4)\] in a paper that employs the optimized Rouse dynamics formalism\]. One usually must resort to approximate methods to include the memory kernel in the analysis. A clear presentation of these methods applied to one-dimensional problems is given in Ref. ppfps:93 .

One way to study the memory function for monomer correlations is through a truncated Mori continued fraction . This becomes rapidly more cumbersome as the order of the expansion is increased. A more convenient expansion of the memory kernel is to first expand the Green’s function in Eq.(28) in the eigenfunctions of $``$, and then expand these eigenfunction in a finite basis (similar to what is done in quantum chemistry calculations) . The basis functions must be chosen with some care for accurate results. Furthermore, in both the continued fraction and eigenfunction expansion methods require equilibrium averages that are difficult to compute for an arbitrary potential, requiring simulations in general. Nevertheless, this is a useful method since the simulation time required for the equilibrium averages can be much shorter than would be required to simulate the time-dependent correlation functions themselves .

To connect to the barrier crossing calculation for protein folding, we need to study the dynamics of the native density, $`\{A\}=\{\rho \}`$. For a general potential, calculating the correlation functions $`\delta \rho _i(t)\delta \rho _j(0)`$ involves the same technical difficulties that arises for the calculation of the monomer correlations. We avoid these complications by approximating the potential by the constrained polymer Hamiltonian, $`H_0`$. Not only is this is a reasonable choice given the interpretation of the order parameters as averages over $`H_0`$, it also allows us to calculate the correlation from Eq.(28) directly.

For an harmonic potential, the propagator in Eq.(28) can be solved analytically . The Gaussian form of $`\rho (𝐫_i)`$ allows the integrals in Eq.(28) to be performed to give the native density correlations (see Appendix for details)

$`\rho _i(t)\rho _j(0)_0`$ $`=`$ $`(detM_{ij}(t))^{3/2}`$ (37)
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \frac{3\alpha ^N}{2a^2}}𝐉_{ij}^TM_{ij}(t)^1𝐉_{ij}\right]`$ (38)
where

$$M_{ij}(t)=\left[\begin{array}{cc}(1+\alpha ^\mathrm{N}G_{ii})& \alpha ^\mathrm{N}G_{ij}(t)\\ \alpha ^\mathrm{N}G_{ij}(t)& (1+\alpha ^\mathrm{N}G_{jj})\end{array}\right],$$
(39)
and, $`𝐉_{ij}`$ is the 2-component difference vector

$$𝐉_{ij}=\left[\begin{array}{c}𝐫_i^N𝐬_i\\ 𝐫_j^N𝐬_j\end{array}\right].$$
(40)
In this expression, we have introduced the time correlation function of monomer positions

$$G_{ij}(t)=\delta 𝐫_i(t)\delta 𝐫_j(0)_0/a^2.$$
(41)
The chain dynamics are contained in $`G_{ij}(t)`$, and the precise form depends on the diffusion matrix $`D_{ij}`$ in Eq.(20). We postpone further discussion of $`G_{ij}(t)`$ until the next section.

The correlations in Eq.(28) for $`\{A\}=\{\rho \}`$ are written as

$$𝒞_{ij}(t)=\delta \rho _i(t)\delta \rho (0)_0=\rho _i(t)\rho _j(0)_0\rho _i\rho _j$$
(42)
where the first term is given by Eq.(37-40) and the second term is found from Eq.(6).

Since approximating the chain potential with $`H_0`$ permits an exact solution for the correlation functions, it may seem as if the discussion of the projected dynamics was an unnecessary detour. However, what we are really trying to calculate is the effective mobility matrix needed to evaluate the folding rate prefactor in Eq.(16).

From Eq.(35), the native density correlations obey the projected relaxation equation

$$_t𝒞(t)=_0^t𝑑t^{}\mu ^\rho (tt^{})\mathrm{\Gamma }^\rho 𝒞(t^{}),$$
(43)
where $`\mathrm{\Gamma }^\rho `$ are the inverse of the static correlations

$$[\mathrm{\Gamma }^\rho ]_{ij}^1=\delta \rho _i\delta \rho _j_0=𝒞_{ij}(0).$$
(44)
Eq.(36) with $`\{A\}=\{\rho \}`$ gives $`\mu ^\rho (t)`$, in terms of projection operators, but we can also deduce an alternative expression for the mobility since $`𝒞(t)`$ is known. Solving Eq.(43) for $`\mu ^\rho (t)`$ by Laplace transforms gives

$$\widehat{\mu }^\rho (\omega )=𝒞(0)\widehat{𝒞}(\omega )^1𝒞(0)\omega 𝒞(0)$$
(45)
To close this expression for the mobility, we only need the Laplace transform of $`𝒞(t)`$ which can be calculated numerically.

As discussed previously, Eq.(43) is consistent with a Gaussian approximation to the free energy

$$\beta F^\rho =\mathrm{const}+\frac{1}{2}\underset{ij}{}\mathrm{\Gamma }_{ij}^\rho \delta \rho _i\delta \rho _j$$
(46)
analogous to the saddle-point expansion of $`\overline{F}[\{\rho ^{}\}]`$ given in Eq.(8). However, even though $`\mathrm{\Gamma }^\rho `$ can be evaluated at the saddle-point constraints in $`\{C^{}\}`$ in $`H_0`$, $`F^\rho `$ is stable against all fluctuations away from $`\delta \rho _i=0`$. This is in contrast to the curvature matrix $`\overline{\mathrm{\Gamma }}^{}`$ that has one unstable direction.

Nevertheless, by approximating the mobility in Eq.(12) by $`\mu (t)\mu ^\rho (t)`$ we can solve the eigenvalue equation Eq.(16) to identify the unstable mode $`𝐮^{}`$ and calculate the prefactor $`\kappa ^{}`$. The last remaining unspecified quantity for the theory is the time dependent monomer correlations $`G(t)`$ which is the subject of the next section.

## IV Chain Dynamics

In this section we specify our model for the chain dynamics, and calculate the monomer pair correlation, $`G(t)`$, with the approximation to the chain potential $`U(𝐑)=H_0`$. This Gaussian approximation simplifies the true microscopic potential that governs conformational changes of the protein chain, and by itself ignores explicit interactions. We attempt to capture some effects of the microscopic dynamics by including frictional forces that are independent of the solvent viscosity, i.e., internal viscosity.

### A Monomer Pair Correlations

Starting from the definition of pair correlations for the Smoluchowski equation \[Eq.(28)\], it is easy to derive an equation for the relaxation of the monomer correlations

$$_t\delta 𝐫_i(t)\delta 𝐫_j(0)=\underset{k}{}\gamma _{ik}^1_kU(\{𝐫(t)\})\delta 𝐫_j(0),$$
(47)
where the friction coefficient is related to the diffusion coefficient by

$$\beta D=\gamma ^1.$$
(48)
In Eq.(47), we have assumed that $`\gamma `$ is independent of $`\{𝐫\}`$, but may still depend on the sequence index. The model for the friction matrix is specified shortly, for now we define a dimensionless friction matrix, $`\stackrel{~}{\gamma }`$:

$$\gamma =\gamma _0\stackrel{~}{\gamma }\gamma _0=6\pi a_{\mathrm{eff}}\eta ,$$
(49)
where we have used Stokes law to define $`\gamma _0`$ with the solvent viscosity $`\eta `$ and typical monomer van der Waals radius $`a_{\mathrm{eff}}`$.

For a given point along the folding route (e.g., a transition state), we assume that the chain dynamics can be described by the reference Hamiltonian. \[Eq.(I-13)\]. Physically, $`H_0`$ corresponds to a polymer confined to the native structure by harmonic constraints with spring constants $`\{3/2a^2C_i\}`$. These constraints act as a non-uniform external field, controlling the fluctuations about the native structure. Taking the gradient of $`H_0`$, and using the definition \[Eq.(I-15)\] of the average monomer position, $`𝐬_i`$, gives

$`{\displaystyle \frac{H_0}{𝐫_i}}`$ $`=`$ $`{\displaystyle \frac{3k_BT}{a^2}}\left[{\displaystyle \underset{j}{}}\mathrm{\Gamma }_{ij}^{(\mathrm{ch})}𝐫_j+C_i(𝐫_i𝐫_i^\mathrm{N})\right]`$ (50)
$`=`$ $`{\displaystyle \frac{3k_BT}{a^2}}{\displaystyle \underset{j}{}}\mathrm{\Gamma }_{ij}^{(0)}\delta 𝐫_j,`$ (51)
where $`\mathrm{\Gamma }_{ij}^{(0)}`$ is given in Eq.(I-17) and $`\delta 𝐫_i=𝐫_i𝐬_i`$ as usual. Substituting this into Eq.(47) gives the matrix equation for the correlations $`G_{ij}(t)=\delta 𝐫_i(t)\delta 𝐫_j(0)/a^2`$,

$$_tG(t)=\sigma \stackrel{~}{\gamma }^1\mathrm{\Gamma }^{(0)}G(t)$$
(52)
where

$$\sigma =\frac{3k_BT}{a^2\gamma _0}=\frac{3D_0}{a^2}.$$
(53)
sets the time scale for the relaxation.

Since Eq.(52) is linear, it can be solved by transforming to normal modes by a similarity transform:

$$Q^1[\stackrel{~}{\gamma }^1\mathrm{\Gamma }^{(0)}]Q=\mathrm{diag}\{\lambda _p\},$$
(54)
where the columns of $`Q`$ are the right eigenvectors of $`[\stackrel{~}{\gamma }^1\mathrm{\Gamma }^{(0)}]`$. Even though $`\stackrel{~}{\gamma }`$ and $`\mathrm{\Gamma }^{(0)}`$ are symmetric, their product is not symmetric in general. Nevertheless, $`Q`$ can diagonalize $`\stackrel{~}{\gamma }`$ and $`\mathrm{\Gamma }^{(0)}`$ separately (though not by a similarity transformation since $`Q`$ is not orthogonal in general) . With

$$Q^T\stackrel{~}{\gamma }Q=\mathrm{diag}\{\nu _p\}Q^T\mathrm{\Gamma }^{(0)}Q=\mathrm{diag}\{\eta _p\}$$
(55)
we have the relationship

$$\lambda _p=\eta _p/\nu _p.$$
(56)
The correlations between monomers, expanded onto the normal modes is

$$G_{ij}(t)=\underset{p}{}Q_{pi}Q_{jp}^T/\eta _p\mathrm{exp}(\sigma \lambda _pt).$$
(57)
We note that equal time correlation functions are the equilibrium correlations $`G_{ij}(t=0)=G_{ij}[\mathrm{\Gamma }^{(0)}]^1`$.

The reference Hamiltonian has two features that distinguish this chain model from the nearest neighbor Rouse chain: chain stiffness, parameterized by $`g`$, and constraints to the native structure, specified by $`\{C\}`$. We consider the effect on the pair correlations of these two contributions separately, to give a point of reference for the model. To simplify this illustration of the dynamics, we consider a (free-draining) diagonal friction matrix, $`\stackrel{~}{\gamma }_{ij}=\delta _{ij}`$. Then, $`\nu _p=1`$ which implies that the relaxation rates are $`\lambda _p=\eta _p`$ and the normal modes are the eigenvectors of the chain connectivity $`\mathrm{\Gamma }^{(\mathrm{ch})}`$.

For an unconstrained polymer ($`C_i=0`$), the dynamics of the stiff, freely rotating chain has been studied quiet extensively . Referring to the definition of $`\mathrm{\Gamma }`$ \[Eq.(I-4)\], $`Q`$ differs from that for the well known rouse modes $`Q^R`$

$$[Q^R]^TK^RQ^R=\mathrm{diag}\{\eta _p^R\}$$
(58)
because of the boundary term, $`g^2/(1g)\mathrm{\Delta }`$, which is only important for very stiff chains . Thus, unless $`g1`$, $`QQ^R`$ for the free-draining, unconstrained polymer:

$$Q_{pi}^\mathrm{R}\mathrm{cos}\left[\frac{\pi p}{n}(i1/2)\right]\eta _p^\mathrm{R}=4\mathrm{sin}^2\left[\frac{\pi p}{2n}\right]$$
(59)
where the mode index takes values $`p=(0,\mathrm{},n1)`$. The relaxation spectrum is then approximately

$$\lambda _p=\eta _p\frac{1g}{1+g}\eta _p^R+\frac{g}{1g^2}(\eta _p^R)^2+B$$
(60)
where the final term is due to the confinement term in Eq.(I-9). The second term of Eq.(60) causes the short wavelength modes (high mode index $`p`$) to relax faster as the chain stiffness increases. In contrast, the relaxation rates for the long wavelength modes (low mode index $`p`$) rates are suppressed by increasing chain stiffness. The number of modes with slowly relaxing rates that are suppressed diminishes as the chain stiffness increases; for the limiting stiffness corresponding to a rigid rod, $`g1`$, the relaxation rates for all modes except one diverge .

When the polymer is constrained by a non-uniform pattern of nonzero constraints $`\{C_i\}`$, the normal modes are not the Rouse modes of a free chain. Fig.1 shows the four slowest modes for a chain with stiffness parameter $`g=0.8`$ (persistence length $`l=5a`$) and various constraints. In addition to the modes for the unconstrained chain shown in Fig.1a, the modes corresponding to selected stationary points along the folding route of $`\lambda `$-repressor are also shown in Fig.1b–f. The slowest modes for the denatured state (Fig.1b) differ only slightly from those for the unconstrained chain (Fig.1a) given by Eq.(59). Referring to the fluctuations shown in Fig.I-9b,the structure at the transition state $`TS_1`$ can be roughly described by saying that residues in helices H4-H5 (indices 53–80) are localized, while the residues in helices H2–H3 (indices 27–47) have the largest fluctuations . These constraints are reflected in the structure of the slowest relaxation modes show in Fig.1c. The motion of the more localized residues is shifted to modes with higher relaxation rates so that these residues are absent in the slowest relaxation modes. In addition, the mode with the slowest relaxation (solid line) is no longer a uniform translation mode, but involves motion of the most unconstrained resides in helices H2–H3 (indices 27–47). The other modes are much the same as those of an unconstrained but shorter chain with one end fixed at the end of helix H4. Similar correspondence between the mode structure and the localization of residues in the transition state ensemble can be seen for the constraints corresponding to $`TS_3`$ and $`TS_4`$ in Fig.1d–e where the modes become increasingly more localized in sequence due to non-sequential constrained residues. The normal modes for the native state (Fig.1f) are more localized with the slowest relaxation involving localized regions with the largest temperature factors.

The mode spectra (relaxation rates) for the free chain, constraints corresponding to transition states, and the native state are shown in Fig.2. As can be seen in this figure, while the spectrum for the short wavelength modes (large mode index) are relatively unaffected, the relaxation rates for all the long wavelength modes (small mode index) become faster as the protein takes on more structure, reflecting the relatively rapid motion toward the native structure as the constraints provided by the native structure become stronger.

### B Internal friction

One can expect that Eq.(52) with harmonic forces can capture the correct physics of the low frequency spectrum of the chain dynamics, but, as in the Debye theory of crystals, the description will not be correct for small length scales which depend on microscopic details. For many polymer physics applications, even though the polymer has a broad range of relaxation rates, the experimentally observed dynamics is dominated by the the long wavelengths (in sequence) motion corresponding to the slowest modes. de Gennes calls the independence of the long wavelength modes on the microscopic conformational changes the “Kuhn Theorem” . Why then is it interesting to incorporate internal viscosity (arising from local activated isomerization) into a model of barrier crossing dynamics? Because we do not know a priori the scale of the motions involved in locally crossing the transition region.

To be more specific, we consider internal friction for the following two reasons: (1) We have seen that the nonuniform constraints to the native structure affects the structure of the normal modes moving them to involve more local motions. It is unclear then how the Kuhn theorem (which is a statement about the delocalized continuum limit) applies here. (2) We are ultimately interested in studying the barrier crossing rate for protein folding, through the Laplace transform of the native density correlations \[Eq.(37)\] which mix different Rouse modes. A correlation function similar to the the native density correlation function arises (for an unconstrained polymer model) in the Wilemski-Fixman theory of diffusion-limited chain closure reactions . In analyzing this theory, Doi showed that for end-to-end monomers, the Laplace transform evaluated at small frequencies is actually dominated by the integral at short times, i.e., all modes contribute to the integral, not just those with the longest relaxation times .

The friction matrix in Eq.(52) $`\gamma =\gamma _0\stackrel{~}{\gamma }`$ is written a sum of two contributions

$$\gamma =\gamma ^\mathrm{S}+\gamma ^{\mathrm{IV}}$$
(61)
where the first term is proportional due to the solvent viscosity, and the second term accounts for the “internal viscosity” and is taken to be independent of the solvent friction. A typical choice for the effective solvent friction is the Oseen tensor from the hydrodynamic mobility . This non-local friction mediated by the solvent damps the motion of one bead due to the motion of another distant beads. It is not difficult to use standard preaveraging approximations to incorporate these hydrodynamic interactions. For simplicity, we neglect this aspect of the problem and consider the free-draining approximation

$$\gamma _{ij}^\mathrm{S}=\gamma _0\delta _{ij}$$
(62)
where $`\gamma _0`$ is defined in Eq.(49).

As previously discussed, some internal friction is due to the neglect of the activated character of the individual conformational transitions of the backbone angles. While there are other ways to incorporate this time scale into a model of chain dynamics , perhaps the most simple correction to the Gaussian description is through the friction matrix. In accordance with the Kuhn theorem, the long wavelength motion should have the usual solvent viscosity dependence being largely unaffected by the microscopic barrier crossing; on the other hand, short wavelengths should tend to relax more slowly, due to the “dynamical stiffness” caused by the slow activation rates for dihedral angle changes. Since the dihedral angle potentials are quite steep these will have a much weaker viscosity dependence.

In the model presented in Ref. erb:mcw:73 , Bazúa and Williams derived an effective frictional force on the $`i^{\mathrm{th}}`$ bead which acts on the relative velocities of adjacent beads:

$$F_{\mathrm{IV}}=\gamma _I(\dot{𝐫}_{i+1}\dot{𝐫}_i)^\mathrm{d}\gamma _I(\dot{𝐫}_i\dot{𝐫}_{i1})^\mathrm{d},$$
(63)
where the superscript denotes the projection of the forces along the bond directions. A common approximation is to write the projected velocities as a sum of the difference between the full velocity vectors $`\{\dot{𝐫}\}`$ and the difference between angular velocity of the beads . The angular component is important for inherently anisotropic phenomena such as flow birefringence, but for our purposes we can neglect the angular rotation component and write the internal friction in the form

$$F_{\mathrm{IV}}=\gamma _IK_{ij}^R_t(\delta 𝐫_j),$$
(64)
where $`K^R`$ denotes the Rouse connectivity matrix \[Eq.(I-5)\], and we have used $`_t(\delta 𝐫_i)=_t𝐫_i`$ since the average monomer position is constant for fixed constraints.

The total friction matrix is the sum of the (free-draining) diagonal solvent friction and the internal friction

$`\gamma _{ij}`$ $`=`$ $`\gamma _0\delta _{ij}+\gamma _IK_{ij}^\mathrm{R}`$ (65)
$`=`$ $`\gamma _0(\delta _{ij}+\gamma _I/\gamma _0K_{ij}^\mathrm{R})`$ (66)
which defines the dimensionless friction matrix $`\stackrel{~}{\gamma }`$ in Eq.(49).

The spectrum of relaxation rates with nonzero internal friction can be illustrated by considering an unconstrained polymer chain ($`C_i=0`$). From Eq.(56) the mode relaxation rates for this case is approximately

$$\lambda _p\frac{\eta _p}{1+\gamma _I/\gamma _0\eta _p^R},$$
(67)
where $`\eta `$ are the eigenvalues of $`\mathrm{\Gamma }^{(0)}`$ (Eq.(60))and $`\eta _p^R`$ are the eigenvalues of the Rouse matrix $`K^R`$. This equation shows that the internal friction leaves the rates of the long wavelength modes (small $`\eta ^R`$) relatively unchanged; this is the Kuhn theorem. The short wavelength modes (larger $`\eta ^R`$), are suppressed by a non-vanishing internal friction. The crossover between the two behaviors, hydrodynamic solvent limited versus conformational interconversion limited, is determined by the value of $`\gamma _I/\gamma _0`$.

For constrained chains, the effect of internal friction on the relaxation modes is different. The internal friction damps rapidly oscillating features of the modes which even occur in the relaxation modes with long relaxation times since they are not identical to Rouse modes. (Note the modes illustrated in Fig.1, which are also different from Rouse mode were calculated for $`\gamma _I=0`$). Nevertheless, for the constraints considered here, this seems to be a minor effect, and the relaxation modes still follow the general behavior of Eq.(67). Fig.3 shows the relaxation spectrum with zero and finite internal friction for the unconstrained chain and for the constraints corresponding to the transition state $`TS_2`$ of the variational calculation for $`\lambda `$-repressor. The spectrum corresponding to the other transition states is similar: the relaxation rates for the slow modes are relatively unaffected, whereas the faster modes have suppressed relaxation rates.

## V Folding Prefactor Example: $`\lambda `$-repressor

The average folding route for the fast folding variant of the $`\lambda `$-repressor protein presented in \[I\] is characterized by sets of constraints $`\{C\}`$ corresponding to local minima and saddle-points of the variational free energy. For each saddle-point, the pair correlations $`G(t)`$ are calculated through the normal mode expansion given by Eq.(57) with the friction matrix given by Eq.(65). These monomer correlations determine the native density correlations $`𝒞(t)`$ through Eq.(42) \[with Eq.(6) and Eq.(37)\], and its Laplace transform $`\widehat{𝒞}(\omega )`$ is taken numerically. With the mobility matrix $`\widehat{\mu }(\omega )`$ given by Eq.(45) and the curvature matrix of the free energy $`\overline{\mathrm{\Gamma }}^{}`$ given by Eq.(10), the eigenvalue equation Eq.(16) is then solved numerically giving the reaction coordinate and folding rate prefactor for that transition state.

In this example, we consider the folding routes obtain for a chain with persistence length of polyalanine, $`l=5a`$ ($`g=0.8`$). The native density $`\{\rho \}`$ also depends on the parameter $`\alpha ^\mathrm{N}`$ controlling the width of the Gaussian measure to the native structure. This width should be small enough to characterize the native structure accurately, but the appropriate value is also limited by the resolution of the model (approximately equal to the root mean square (RMS) fluctuations of the bond length $`a`$). For large $`\alpha ^\mathrm{N}`$ (small width), $`\{\rho \}`$ has low values even for the native state since the average positions are not precisely equal to the native structure at finite temperature; for $`\alpha ^\mathrm{N}=1`$ the average value of $`\rho _i`$ evaluated at the native state is 0.3. For the native densities shown in Fig.4 with $`\alpha ^\mathrm{N}=0.5`$, the average of $`\rho _i`$ at the native state increases to $`\rho _i`$ is 0.5. We consider this to be a reasonable value to use for $`\alpha ^\mathrm{N}`$. For smaller $`\alpha ^\mathrm{N}`$ (larger widths), the separation between $`\{\rho \}`$ evaluated at the native and globule state is more distinct, though unfolded regions of the transition states have native densities closer to the native values. The dependence of the prefactor on the width $`\alpha ^\mathrm{N}`$ is considered below.

The local reaction coordinates (unstable modes) for the four transition states along the folding route are shown in Fig.5. These unstable modes are quite understandable in terms of the structure inferred by the temperature factors in Fig.I-9b.The unstable mode of $`TS_1`$ (Fig.5a) is rather delocalized, consisting of residues in helices H4–H5 and modulates with a wavelength approximately equal to the persistence length $`l/a=5`$; the partial localization of helix H1 indicated by the temperature factors is not a part of the normal mode. The mode for $`TS_2`$ (Fig.5b) has a significant component of formation of helices H4–H5 as in $`TS_1`$, but the localized region around a residue in H1 is a much larger component. The unstable mode at $`TS_3`$ is also localized, but consists mainly of a few residues separated in sequence in helix H1 and near helix H4. Finally, the unstable mode at $`TS_4`$ is delocalized consisting of residues in helices H2–H3 which are the last helices to form. The features of the reaction coordinates resemble the unstable modes of the free energy ($`\overline{\mathrm{\Gamma }}`$), though the relative amplitudes of the peaks in a given reaction coordinate are altered by the chain dynamics through the mobility $`\widehat{\mu }(\omega )`$.

Since the local reaction coordinates along the average folding route are very different, one from each other, a reduced description of the kinetics based on a single reaction coordinate can not be used when dynamics is examined at this level of detail. Attempts to find perfect global dynamical coordinates in lattice simulations of perfectly funneled models are consistent with this . In the present description, each dynamical reaction coordinate applies in the vicinity of one saddle-point. The changes of the local reaction coordinate along the folding route are intimately related to the fine structure of the free energy profile. Notice that the effects on predicted reaction rates are minor, however. Surmounting the highest barrier gives the rate within a factor of two (vide supra). Calculations of one-dimensional barrier crossing with shallow local minima (intermediates) also show why having a perfectly good local detailed description of the free energy profile need not yield a strictly defined global reaction coordinate . At the same time, a much cruder reaction coordinate can function perfectly well for predicting rates if the appropriate configurational diffusion coefficient is used .

The prefactor $`k_0=|\kappa ^{}|/2\pi `$ for each transition state is plotted as a function of the internal friction $`\gamma _I/\gamma _0`$ in Fig.6. When the unstable curvature of the free energy is large, there is less re-crossing at the top of the barrier leading to a larger prefactor. For small $`\gamma _I/\gamma _0`$ the prefactors increase in the same order as the corresponding curvatures of the free energy ($`\overline{\mathrm{\Gamma }}_0`$): $`k_0^{(4)}<k_0^{(1)}<k_0^{(2)}<k_0^{(3)}`$. The relative values of the prefactors at different transition states depend on the chain dynamics through the mobility matrix $`\widehat{\mu }(\omega )`$, however, and do not follow the proportions of the corresponding free energy curvatures. For small values of $`\gamma _I/\gamma _0`$, the prefactor is proportional to the bare diffusion coefficient $`D_0=k_BT/\gamma _0`$. The prefactors decrease with increasing $`\gamma _I/\gamma _0`$ because of the suppressed relaxation rates of the fastest modes. This turnover occurs around $`\gamma _I/\gamma _0=O(1)`$, though it is different for each transition state.

To put these results in laboratory units for polypeptides, we need to specify the monomer diffusion coefficient $`D_0`$ and the internal viscosity parameter $`\gamma _I/\gamma _0`$. We estimate the typical hydrodynamic radius in Eq.(49) to be the sum of the monomer spacing $`a=3.8`$Å and the average van der Waals radius of a typical side chain $`a_{\mathrm{vdw}}3`$ Å : $`a_{\mathrm{eff}}7`$Å. Using the viscosity of water at room temperature $`\eta _0=1`$ centipoise gives the a monomer diffusion coefficient $`D_0^{\mathrm{H}_2\mathrm{O}}=3\times 10^6\mathrm{cm}^2/\mathrm{s}`$. We estimate the internal viscosity by the elementary bond flipping rate of the polymer backbone. Calculations of the $`\beta \alpha `$ conformations of di-alanine give a value for $`\tau _{\mathrm{flip}}^110\mathrm{n}\mathrm{s}^1`$ . One can expect this rate to be reduced when the flipping transition occurs in a long polypeptide chain since the neighboring monomers must adjust to some extent in order to minimize the motion of the chain through the solvent . Judging from studies on model chains one expects about a factor of three slower, we estimate a typical value to be citejs:eh:80,eh:js:82 $`\tau _{\mathrm{flip}}^1=3\mathrm{ns}^1`$ . It would be nice for simulators to pin this number down for the peptide isomerization rate in long chains. It is obvious that the relaxation rates describing the chain dynamics should not exceed this elementary rate. Associating the internal friction with this isomerization rate, we set $`\gamma _I/\gamma _0`$ so that the free chain relaxation rates of the fastest mode equals $`\tau _{\mathrm{flip}}^1`$. For a chain with persistence length $`l=5a`$ and monomer diffusion coefficient $`D_0^{H_2O}`$, we find $`\gamma _I/\gamma _0=9.1`$.

In Fig.7, the prefactor $`k_0`$ is plotted against the inverse monomer diffusion coefficient $`D_0^1`$ in laboratory units. The internal friction is set to the fixed value $`\gamma _I=9.1\gamma _0`$ for $`\gamma _0/k_BT_\mathrm{f}=(D_0^{\mathrm{H}_20})^1`$. The prefactor is proportional to the solvent viscosity ($`D_0^1\eta _0`$) for large viscosities and saturates at low viscosities. This turnover is reminiscent of the viscosity dependence of the prefactor from Kramer’s theory , but here the turnover is driven indirectly by the internal friction of the chain rather than explicit inertial terms reflecting an energy limited regime. The prefactor in pure water, indicated by the dashed line in Fig.7, is still in the linear regime, but near the turnover. The prefactors of the different transition states for $`D_0=D_0^{\mathrm{H}_2\mathrm{O}}`$ (given in Table I) are on the order of $`10\mu \mathrm{s}^1`$. As a point of reference, the unstable mode frequency for $`TS_2`$, $`|\kappa ^{}|=6.9\times 10^3\sigma `$, is approximately equal to the relaxation rate of the $`p=6`$ normal mode of a free chain of $`n=80`$ monomers and persistence length $`l=5a`$ ($`\lambda _{p=5}=7.2\times 10^3`$). The normal mode for this motion has a node to node period of 20 monomers, twice the length used in the estimates based on the smallest loop closure .

The folding rate includes the activation barrier (calculated with this model in \[I\]) as well as the prefactor. The calculated thermodynamic and kinetic quantities needed for the folding rate are given in Table II for two thermodynamic conditions. Adjusting the energy scale of the contacts $`ϵ_0/k_BT`$ fixes the equilibrium constant $`K=k_f/k_u`$ (the ratio of the folding and unfolding rates), so that $`\mathrm{ln}K`$ gives the difference in free energy between the globule and native state (in units of $`k_BT`$). For the very stabilizing condition $`\mathrm{ln}K=7.8`$, the prefactor is $`k_0=2.4\mu \mathrm{s}^1`$ and the barrier is $`\mathrm{\Delta }F^{}/k_BT=3.4`$, giving a folding time of $`1/k=13\mu \mathrm{s}`$. This agrees well with the experimental rate $`1/k=12\mu \mathrm{s}`$ at this stability (extrapolated to 0M denaturant) . At the transition midpoint $`\mathrm{ln}K=0`$, the calculated prefactor is $`7.6\mu \mathrm{s}`$, and with the barrier height $`\mathrm{\Delta }F^{}/k_BT=5.1`$, we get a folding time of $`1/k=23\mu \mathrm{s}`$ which is considerably faster than the measured folding time of $`1/k10m\mathrm{s}`$ . This suggest the microscopic calculation for the perfectly funneled surface under estimates the barrier by about $`5k_BT`$ at the folding transition temperature.

The calculated inverse prefactor for $`\lambda `$-repressor at the transition midpoint $`1/k_00.1\mu \mathrm{s}`$ and under highly stabilizing condition $`1/k_00.5\mu \mathrm{s}`$ are both faster than the estimated fastest possible folding time  $`\tau 1\mu \mathrm{s}`$. The latter is the time it takes to form a typical contact as determined by the fluorescence or chemical quenching rate of a pair of residues separated by a typical loop length. It is very difficult to imagine that folding can occur on times faster than contact formation, and the introduction of this speed limit is important conceptually as a reference timescale for protein folding dynamics. It is tempting to identify the speed limit with the folding prefactor since the fastest rate is obtained from Eq.(1) for a vanishing free energy barrier. But this is not precisely correct since the separation of time scales assumed in activated kinetics is not true for small barriers. As seen in Table II, reducing the activation barrier (by stabilizing the native state, for example) changes the prefactor since the saddle–point curvature of the free energy changes. If the prefactor at both stabilities in this example were $`1/k_0=0.1\mu \mathrm{s}`$, the folding time with the reduced activation barrier would be a factor of approximately 6 times smaller than at the transition midpoint, or $`1/k4\mu \mathrm{s}`$, which approaches the speed limit. The prefactor is found to be approximately three times smaller than at the midpoint, however, resulting in a larger folding time $`1/k12\mu \mathrm{s}`$. Pushing the conditions to the limit of a vanishing barrier, it could happen that the barrier is very low but still has a large unstable mode free energy curvature, resulting in a calculated folding time somewhat faster than the speed limit. For such near downhill folding , however, the folding would not be single exponential process dominated by barrier crossing, but involves kinetics more directly connected to the time for contact formation than the time reflected in the apparent prefactor.

It is interesting to consider the dependence of the prefactor on the chain dynamical parameters since the values for the monomer diffusion coefficient and the isomerization rate in a polymer chain are not known precisely (see next section). Conservatively, the values we used are probably accurate to within a factor of two or three. As long as the internal friction is not too large, the monomer diffusion coefficient is the more important parameter for the absolute prefactor (in $`\mathrm{s}^1`$); since the time scale for the monomer relaxation is $`3D_0/a^2`$, the absolute prefactor is proportional to this scaling. Table III gives the prefactor from $`TS_2`$ using the free alanine diffusion coefficient and the di-alanine isomerization rate. As shown in Table III, the prefactor varies by about a factor of three when $`D_0^{\mathrm{H}_2\mathrm{O}}`$ is the diffusion coefficient of a free alanine molecule. This is the same factor as the ratio of the monomer diffusion coefficient, though it also depends on the value of the internal friction (i.e., $`\tau _{\mathrm{flip}}^1`$) as well.

The value of the Gaussian width $`\alpha ^\mathrm{N}`$ defining the native density is a parameter of the model, not a measurable quantity, so we must also consider the sensitivity of our results to the value of $`\alpha ^\mathrm{N}`$. As shown in Fig.8, the prefactors $`k_0`$ are approximately linear in $`\alpha ^\mathrm{N}`$, increasing by a factor of about two for $`TS_3`$ and less than a factor of unity for the other transition states as $`\alpha ^\mathrm{N}`$ varies from 0.25 to 1.0. The different slopes are reflect the different increases in curvatures of the free energies as the Gaussian width narrows. This degree of uncertainty seems reasonable, but the different slopes suggests that the dependence is a delicate issue, depending on the detailed structure of each transition state.

## VI Concluding Remarks

In this paper, we derived a microscopic theory for the dynamics of protein folding barrier crossing when non-native contacts and trapping effects can be neglected. The variational theory presented in \[I\] characterizes the average folding routes as connected minima and transition states. These routes are specified through variational parameters constraining the protein chain to inhomogeneously order about the native structure. The quadratic reference Hamiltonian from this theory along with forces from solvent and internal friction govern the chain dynamics near the saddle-points. The barrier crossing dynamics is calculated from a generalization of both Grote-Hynes and Langer rate theories to include a frequency dependent friction in multi-dimensions. The memory friction itself is approximated by linear response theory and exploits the harmonic potential of the model. The theory gives local reaction coordinates and folding rate prefactors for specific proteins with known native structure.

For the calculated average folding route of $`\lambda `$-repressor, we find a folding prefactor $`k_010\mu \mathrm{s}^1`$ at the transition midpoint and $`k_01\mu \mathrm{s}^1`$ under highly stabilizing conditions. Using the free energy barriers from this model, the calculated folding rate agrees well with the measured rate extrapolated to zero denaturant, but not at the transition midpoint. Considering other reasonable values for the parameters in the model suggests that the calculated prefactor may be off by a factor of about three. These calculations also suggest that the folding rate prefactor near water solvent conditions scales linearly with solvent viscosity. This dependence agrees with experimental studies on other two state proteins  and simulations of off-lattice model proteins . Nevertheless, these conditions are near a turnover regime for which the internal friction of the polymer chain dominates, ultimately causing the prefactor to saturate upon decreasing solvent viscosity. Thus folding of polymers with other backbones may be quite different.

The discrepancy between the calculated and measured folding rate at the transition midpoint $`(T_\mathrm{f})`$ suggests that the model underestimates the activation barrier or overestimates prefactor for this stability. Since the model has many simplifications that could influence both the thermodynamics and the dynamics, it is impossible without further investigation to identify which primarily needs improvement. It may be combination of both.

One aspect of the present model is that the native state is quite a bit more flexible than observed in the laboratory, with RMS fluctuations on the order of the distance between adjacent monomers $`a`$ (the spatial resolution of the the harmonic chain). In comparison, temperature factors from x-ray crystallography correspond to RMS fluctuation of less than 1Å $`(a/4)`$. A more realistic non-Gaussian polymer backbone would increase the resolution of the model. Such a non-Gaussian model could still be treated using the variational formalism. Still retaining the Gaussian description of the chain, multi-body forces arising from neglected degrees of freedom would also tend to increase the rigidity of the native state and are a likely contributor. Recent off-lattice simulations have shown that effective multi-body forces sharpen interfaces between folded and unfolded regions and increase activation barriers . As dicussed in \[I\], contiguous sequence approximations  also give inherently sharp interfaces and larger activation barriers even for pairwise additive forces since the conformations of the residues are restricted to be either folded or unfolded. (The contiguous sequence approximations with multi-body forces would give even higher barriers.) If multi-body forces are indeed important, the completely folded or unfolded conformations assumed in these approximations may well be more appropriate than expected from simulations of lattice and off-lattice models with only pairwise forces. Still, the interfaces are probably not as cleanly sharp as envisioned in the contigous sequence approximations. Using the calculated prefactor at the transition midpoint, $`k_08\mu \mathrm{s}^1`$, and the measured folding rate of $`100\mathrm{s}^1`$ gives a barrier of approximately $`11k_BT_\mathrm{f}`$. This value is greater than both the activation barrier from the variational theory $`(\mathrm{\Delta }F^{}5k_BT_\mathrm{f})`$ and the double sequence approximation calculated in \[I\] $`(\mathrm{\Delta }F^{}8k_BT_\mathrm{f})`$. Multi-body forces may sharpen the distinction between folded and unfolded regions in the transition state ensemble to give an activation barrier consistent with the prefactor calculation at the transition midpoint.

Alternatively, it is possible the discrepancy at $`T_\mathrm{f}`$ arises because of neglected dynamical effects that could slow the chain dynamics, reducing the folding rate prefactor. In this paper, the internal friction is interpreted as a correction to account for the activated isomerization reaction of the backbone. As mentioned in the introduction, there are other possible contributors to this friction. For example, non-native interactions are presently neglected in the model. These interactions contribute to the memory friction and therefore the internal viscosity. They can be treated through mode coupling calculations similar to those we have done for a random heteropolymer  (see also Refs. jrr:eis:96 ). A more careful treatment of the solvent may also be necessary for an accurate description of the chain dynamics. Detailed calculations of the reaction dynamics for isomerization transition in di-alanine (a much less complicated system) suggest the solvent may be an essential component of even the local reaction coordinate . Furthermore, contact formation dynamics may also be slowed because contact formation excludes the solvent, thereby adding an additional activation barrier and timescale into the dynamics. Supporting this possibility, we note that simulations of minimal protein folding models with desolvation barriers in the pair potential bring simulations into qualitative agreement with pressure studies on folding rates .

Since we expect some of these neglected internal friction contributions to influence protein folding dynamics, we have been careful in this paper to use “bare” parameters (and investigate other reasonable estimates) instead of treating the them as fitting parameters in the theory. This is an attempt to avoid missing a signal that such new dynamical considerations are indeed important by using inappropriate effective parameters. The folding rate alone, however, is not enough to determine the accuracy of the calculated prefactor, since the barrier height may not yet be accurately enough detetermined.

Consequently, these proposed effects are perhaps best studied at this point using more direct experimental probes such as intrachain fluorescence quenching measurements. Fluorescence quenching is a powerful tool to study the average structure of ensembles  and has recently been used as a probe for experiments on single molecules that allow distributions to be addressed . Motivated primarily by the relevance to protein folding dynamics, some fluorescence quenching experiments have focused on less complicated systems such as short peptides . These measurements of simple polymer pair dynamics can guide the development of improved theoretical models helpful to understanding protein folding kinetics. This is particularly evident in the present model, since the expression for the effective mobility \[Eq.(45)\] is formally related to approximations to the intrachain quenching rate in polymers .

Experiments monitoring the intrachain contact formation through fluorescence quenching suggest that the dynamical issues mentioned above play a role in the pair contact dynamics as well. Fitting data to effective models for the pair dynamics, the apparent diffusion coefficient comes out to be 10–20 times smaller than one would expect from the diffusion coefficients of the free probes . In those models, the polymer is treated as a single pair connected by an entropic spring controlling the mean square fluctuation of chain ends  where the effective diffusion coefficient of the pair, $`D_{\mathrm{eff}}=2D_0`$, is a short time approximation to the full dynamics including all the monomers . Haas has suggested that the internal friction of the chain may explain the descrepancy between $`D_{\mathrm{eff}}`$ and the fitted value , but no detailed calculation has been done. This is a problem we wish to return to . Doi’s arguments suggest that the short time properties may dominate . We note that the inclusion of internal viscosity from the bond flipping rate as presented here $`(\gamma _I/\gamma _010)`$ reduces the effective short time diffusion coefficient to a value $`D_{\mathrm{eff}}D_0/2`$, narrowing the discrepancy to within a factor of about 2–5. Since the validity of this short time approximation can be questioned (e.g., the rate depends on the quenching volume ), we feel a solution resolution to this puzzle requires a more complete calculation.

The primary experimental technique used to study the structure of the transition state ensemble is the $`\varphi `$-value analysis developed by Fersht . By comparing the change in the folding rate to the change in stability upon mutation, one can infer the extent to which that site is structured in the transition state ensemble. To interpret these experiments, one assumes that the mutation alters the kinetics through a change in the free energy profile, leaving the prefactor unchanged. Sometimes, however, the kinetics and thermodynamics are seen not to be directly coupled. This may well be a sign of frustration from non-stabilizing native contacts or non-native contacts influencing the kinetics of the transition . It is also possible however that these mutations affect the chain dynamics by modifying the internal friction at that residue. It is easy to extend the present calculation to indicate how the internal friction could be exploited to find “dynamical mutants” in which a residue is replaced by one with higher dihedral barriers (i.e., a residue with larger local internal friction) perturbing the prefactor, without altering the stability. Proline isomerization coupled to folding are an example of this, but other nonnatural amino acids may be still more promising probes of local reaction coordinates . Mutations of the n-src loop of the SH3 domain have very small effect on the equilibrium constant (preventing the $`\varphi `$-value analysis), but still impact the folding rate and thus may be examples of dynamical mutants . More extreme alterations to the sequence have been necessary such as lengthening the loop length in order to characterize this region of the transition state ensemble . A better understanding of the chain dynamics should help suggest additional and more direct probes of similar situations using dynamical mutants.

## Acknowledgments

This work has been supported by NIH Grant No. PHS 2 R01 GM44557.

## A

In this appendix, we outline the calculation of the native density correlations $`\rho _i(t)\rho _j(0)_0`$ given by Eq.(37).

Since $`H_0`$ is harmonic, the Green’s function for the diffusion operator is known and has the form of a multi-variable Gaussian containing all the pair correlations . The native density correlations only depend on two monomer positions for a fixed $`i`$ and $`j`$, so that Eq.(28) can be simplified by introducing delta functions $`𝑑𝐫_1𝑑𝐫_2\delta (𝐫_1𝐫_i)\delta (𝐫_2𝐫_j)`$ and integrating over all monomer variables $`\{𝐫\}`$. This reduces Eq.(28) to a double integral over $`𝐫_1`$ and $`𝐫_2`$, and the reduced Green’s function contains only the monomer correlations of the monomer pair, as can be expected for a Gaussian process.

It is perhaps easier to make use this Gaussian property from the beginning, and obtain the same result as using the full Green’s function as described above.

Consider a fixed $`i`$ and $`j`$ pair. First we define some notation. We define the following 2-component vectors:

$$𝐑=\left[\begin{array}{c}𝐫_i(t)\\ 𝐫_j(t^{})\end{array}\right]𝐒=\left[\begin{array}{c}𝐬_i\\ 𝐬_j\end{array}\right]𝐑^\mathrm{N}=\left[\begin{array}{c}𝐫_i^\mathrm{N}\\ 𝐫_j^\mathrm{N}\end{array}\right]$$
(A1)
where, $`t`$ and $`t^{}`$ denote two times, and e.g., $`𝐬_i=r_i_0`$ and $`𝐫_i^\mathrm{N}`$ is the average and native position vector of the $`i^{\mathrm{th}}`$ monomer, respectively. Let $`\chi (tt^{})^1`$ be the symmetric $`2\times 2`$ matrix of monomer correlations $`\chi (tt^{})^1=\delta 𝐑(t)\delta 𝐑(t^{})^T_0/a^2`$, with $`\delta 𝐑(t)=𝐑(t)𝐒`$:

$$\chi (|tt^{}|)^1=\left[\begin{array}{cc}G_{ii}& G_{ij}(|tt^{}|)\\ G_{ij}(|tt^{}|)& G_{jj}\end{array}\right]$$
(A2)
To construct $`\chi (|tt^{}|)`$, we have denoted the equal time correlation function by the static correlations $`G_{ii}`$ and $`G_{jj}`$, and used symmetry of the correlations $`G_{ij}(t)=G_{ji}(t)`$.

Since this is a Gaussian process, the correlation function $`\rho _i(t)\rho _j(t)_0`$ can be written as

$$\rho _i(t)\rho _j(t^{})_0=_{𝐫_i,𝐫_j}\rho _i(𝐫_i)\rho _j(𝐫_j)\stackrel{~}{P}(\delta 𝐫_i,\delta 𝐫_j;|tt^{}|),$$
(A3)
where

$`\stackrel{~}{P}(\delta 𝐫_i,\delta 𝐫_j;|tt^{}|)={\displaystyle \frac{3}{2\pi a^2}}(det\chi (|tt^{}|))^{3/2}`$ (A5)
$`\times \mathrm{exp}\left[{\displaystyle \frac{3}{2a^2}}\delta 𝐑^T\chi (|tt^{}|)\delta 𝐑\right].`$
In the rest of the derivation, we suppress the time dependence to simplify notation. Substituting Eq.(4) for $`\rho (𝐫_i)\rho (𝐫_j)`$ and shifting the integration variables to $`\delta 𝐑`$ ($`𝐑=\delta 𝐑+𝐒`$) leads to

$`\rho _i(t)\rho _j(t)_0={\displaystyle \frac{3}{2\pi a^2}}(det\chi )^{3/2}\mathrm{exp}\left[{\displaystyle \frac{3\alpha ^\mathrm{N}}{2a^2}}(𝐑^\mathrm{N}𝐒)^2\right]`$ (A7)
$`\times {\displaystyle _{\delta 𝐫_i,\delta 𝐫_j}}\mathrm{exp}[{\displaystyle \frac{3}{2a^2}}\delta 𝐑\chi ^{}\delta 𝐑+{\displaystyle \frac{3}{a^2}}𝐉\delta 𝐑],`$
where

$$\chi ^{}=\chi +\alpha 𝐈I=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$$
(A8)
$$𝐉=𝐑^\mathrm{N}𝐒.$$
(A9)
Performing this integral by completing the square in the exponent gives

$`\rho _i(t)\rho _j(t)_0`$ $`=`$ $`\left[{\displaystyle \frac{det\chi }{det\chi ^{}}}\right]^{3/2}\mathrm{exp}\left[{\displaystyle \frac{3\alpha ^\mathrm{N}}{2a^2}}𝐉^2\right]`$ (A10)
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \frac{3(\alpha ^\mathrm{N})^2}{2a^2}}𝐉^T\chi ^1𝐉\right].`$ (A11)
To simplify this equation, it is convenient to re-write $`\chi ^{}`$ as

$$\chi ^{}=\chi M$$
(A12)
where

$`M`$ $`=`$ $`(I+\alpha ^\mathrm{N}\chi ^1)`$ (A13)
$`=`$ $`\left[\begin{array}{cc}1+\alpha ^\mathrm{N}G_{ii}& \alpha ^\mathrm{N}G_{ij}(|tt^{}|)\\ \alpha ^\mathrm{N}G_{ij}(|tt^{}|)& 1+\alpha ^\mathrm{N}G_{ij}\end{array}\right].`$ (A16)
Then we can express $`\chi ^1`$ as

$$\chi ^1=M^1\chi ^1.$$
(A17)
Now since

$`I`$ $`=`$ $`M^1M`$ (A18)
$`=`$ $`M^1(I+\alpha ^\mathrm{N}\chi ^1)`$ (A19)
$`=`$ $`M^1+\alpha ^\mathrm{N}M^1\chi ^1`$ (A20)
$`=`$ $`M^1+\alpha ^\mathrm{N}\chi ^1,`$ (A21)
we can express $`\chi ^1`$ as

$$\chi ^1=\frac{1}{\alpha ^\mathrm{N}}(IM^1).$$
(A22)
Consequently, the combined argument the exponential factors in Eq.(A10) is

$$\frac{3(\alpha ^\mathrm{N})^2}{2a^2}(𝐉^T\chi ^1𝐉\frac{1}{\alpha ^\mathrm{N}}𝐉^2)=\frac{3\alpha ^\mathrm{N}}{2a^2}𝐉^TM^1𝐉.$$
(A23)
The determinate factor can be simplified using

$$det\chi ^{}=det\chi detM.$$
(A24)
Thus, Eq.(A10) becomes

$$\rho _i(t)\rho _j(t)_0=(detM)^{3/2}\mathrm{exp}\left[\frac{3\alpha ^\mathrm{N}}{2a^2}𝐉^TM^1𝐉\right],$$
(A25)
which is Eq.(37).