text:
# Elementary excitations in one-dimensional spin-orbital models: neutral and charged solitons and their bound states

## I Introduction

One-dimensional systems with coupled spin and orbital degrees of freedom have been attracting a considerable attention during the last two years. This interest is partly motivated by the progress in the experimental study of the quasi-one-dimensional spin-gap materials $`\mathrm{Na}_2\mathrm{Ti}_2\mathrm{Sb}_2\mathrm{O}`$ and $`\mathrm{NaV}_2\mathrm{O}_5`$, which are believed to be described by a two-band orbitally degenerate Hubbard model at quarter filling. Orbital degrees of freedom may be described with the help of pseudospin-1/2 variables; the corresponding effective Hamiltonian for the two-band Hubbard model was derived long ago by Kugel and Khomskii  and has the following form:

$`\widehat{H}`$ $`=`$ $`K{\displaystyle \underset{n}{}}(𝑺_n𝑺_{n+1})(𝝉_n𝝉_{n+1}+\epsilon ^{}\tau _n^z\tau _{n+1}^z)`$ (1)
$`+`$ $`{\displaystyle \underset{n}{}}\left\{J_s(𝑺_n𝑺_{n+1})+J_\tau (𝝉_n𝝉_{n+1}+\epsilon \tau _n^z\tau _{n+1}^z)\right\},`$ (2)
where $`𝑺_n`$ and $`𝝉_n`$ are respectively the spin-$`\frac{1}{2}`$ and pseudospin-$`\frac{1}{2}`$ operators at the site $`n`$. The Hamiltonian (1) can be viewed as two spin-$`\frac{1}{2}`$ chains coupled by a biquadratic interaction only (a generalized spin ladder).

Generally, the model has an $`\mathrm{SU}(2)\times \mathrm{U}(1)`$ symmetry, so that the total spin $`(S,S^z)`$ of the $`𝑺`$-chain, as well as the $`z`$-projection $`\tau ^z`$ of the total spin of the $`𝝉`$-chain are good quantum numbers. Under certain simplifying assumptions (neglecting Hund’s rule coupling, nearest neighbor hopping between the same type of orbitals only, and only one Coulomb on-site repulsion constant), one obtains the above Hamiltonian with $`\epsilon =\epsilon ^{}=0`$, $`J_s=J_\tau =\frac{1}{4}K`$ which possesses the hidden SU(4) symmetry.  At this special point, the model is Bethe ansatz solvable  and gapless. Depending on the microscopic details of the interaction, this high symmetry can be broken in several ways. For example, finite Hund’s rule coupling or existence of more than one Coulomb repulsion constant makes $`J_s`$ and $`J_\tau `$ different, and local crystal fields can induce considerable anisotropy in the orbital sector.

The phase diagram of the isotropic ($`\epsilon =\epsilon ^{}=0`$) $`\mathrm{SU}(2)\times \mathrm{SU}(2)`$ version of the model (1) in the vicinity of the $`SU(4)`$ point was analytically studied recently by several authors.  Extensive density matrix renormalization group (DMRG) studies of the phase diagram  established the form of phase boundaries and the quantum numbers of the lowest-lying excitations.

Moving away off the $`SU(4)`$ point towards larger $`J_s`$, $`J_\tau `$, one runs into the spontaneously dimerized phase with a finite gap and twofold degenerate ground state. The structure of the elementary excitation spectrum in the dimerized phase is known rather incompletely. The weak coupling ($`J_{s,\tau }K`$) region of the dimerized phase is known to be a realization of the so-called non-Haldane spin liquid  where magnons become incoherent excitations since they are unstable against the decay into soliton-antisoliton pairs. Variational study for a special “dimer point” $`J_s=J_\tau =\frac{3}{4}K`$, where the exact ground state in a form of a simple valence bond solid state is known,  confirms that solitons are the lowest excited states also in case of strong coupling. However, recent DMRG results  revealed a complicated pattern of the change in quantum numbers of the lowest-lying excited states in the vicinity of this special point; the standard DMRG technique provides only the information about the quantum numbers and energy of the lowest excited state, but not the information on the entire spectrum. Recent work  using the so-called dynamical DMRG method contains data for the dynamical structure factor, but for only one point in the dimerized phase, namely for the dimer point.

Further, realistic material-relevant models are generically anisotropic ($`\epsilon ,\epsilon ^{}0`$) and asymmetric ($`J_sJ_\tau `$);  the anisotropic model of the same type arises also as an effective Hamiltonian for a spin tube. Thus we think it is of interest to get insight into the properties of the elementary excitations in the dimerized phase of the more general asymmetric/anisotropic model, which is the aim of the present paper.

In section II we introduce briefly the “dimer line” in the space of model parameters where the exact twofold degenerate spontaneously dimerized ground state is known. Based on this picture of the ground state, in section III we present the variational study of the interplay between different types of elementary excitations, namely, between soliton-antisoliton pairs and their bound states (magnons). We show that soliton pairs are the lowest excitations only in a narrow region of the phase diagram, and away from the dimer line the low-energy physics is determined by magnon-type localized excitations which are essentially soliton-antisoliton bound states. An additional phase boundary is shown to appear when the anisotropy exceeds a certain critical value. Variational dispersion curves are found to agree qualitatively with the exact diagonalization (Lanczos) results in the vicinity of the dimer line.

In section IV we describe some exact results for the doped spin-orbital model on the “dimer line”. We show that the natural charge excitation in this case is a “charged soliton”, i.e., the bound state of a hole and a soliton in the dimer order. It is found that under certain condition “charged solitons” themselves exhibit a tendency to form bound states; exact expressions for the energy of those bound states are presented.

## II Exact ground states

In addition to the $`SU(4)`$ model, there are several other particular points in the phase diagram of (1), for which exact results are available. Particularly, at the point $`J_s=J_\tau =3K/4`$ the exact ground state  is a product of checkerboard-ordered spin and orbital singlets (see Fig. 1), so that the spin and orbital degrees of freedom are completely decoupled in the ground state. This ground state is spontaneously dimerized and thus twofold degenerate in case of periodic boundary conditions (for open boundaries the degeneracy would be 8-fold).

One can show  that the same ground state persists in a finite range of anisotropy in the case of the special choice $`\epsilon =\epsilon ^{}`$. In this case, up to a constant term, the model (1) can be recast in the form

$$\widehat{H}=\underset{n}{}(𝑺_n𝑺_{n+1}+J_\tau )(𝝉_n𝝉_{n+1}+\epsilon \tau _n^z\tau _{n+1}^z+J_s),$$
(3)
where we have set $`K`$ to be the unit energy scale. Then it is easy to see that the choice $`J_\tau =3/4`$, $`J_s=(3+\epsilon )/4`$ makes the checkerboard-type singlet product wave function an eigenstate, and one can prove that it stays the ground state as long as $`2<\epsilon <\mathrm{}`$. Thus, in the 3D phase space $`(J_s,J_\tau ,\epsilon )`$ one has a line with an exactly known ground state; in what follows we will call it the “dimer line” for the sake of brevity. This line ends at the multicritical point $`\epsilon =2`$, where the ground state becomes degenerate with infinitely many other states: e.g., changing any of the orbital ($`𝝉`$) singlets into a triplet with $`\tau ^z=\pm 1`$ does not change the energy. The point $`\epsilon =+\mathrm{}`$ is another multicritical point where the triplets with $`\tau ^z=0`$ can be created with no energy cost. For $`\epsilon <2`$ the ground state is also exactly known  and is a product of the ferromagnetic state in $`𝝉`$-sector and the antiferromagnetic Bethe-ansatz state in the $`𝑺`$-sector. The energy of ferromagnons softens at $`k=\pi `$ when $`\epsilon `$ tends to $`2`$.

We have also checked that (3) at $`J_\tau =3/4`$, $`J_s=(3+\epsilon )/4`$ is the only nontrivial model with exactly known ground state (within the physically relevant subspace of Hamiltonians defined by (1)) which can be found within the matrix-product based optimum ground states approach (see e.g. Ref.  and references therein), at least if one uses a rather general matrix-product ansatz

$`|\mathrm{\Psi }=\text{Tr}(g_1(C)g_2(\stackrel{~}{C})\mathrm{}g_{N1}(C)g_N(\stackrel{~}{C}))`$

where $`C\{x,y,x^{},y^{},p,q\}`$, and the elementary matrix $`g_n(C)`$ has the following form:

$`g_n(C)=\left(\begin{array}{cc}x|_n+y|_n\hfill & \hfill p|_n\\ q|_n\hfill & \hfill x^{}|_n+y^{}|_n\end{array}\right).`$

Here $`n`$ denotes the rung of the “ladder,” and the first and second arrows denote the states of $`𝑺`$ and $`𝝉`$ spins, respectively.

## III Variational study of the elementary excitations

In the vicinity of the dimer line, the elementary excitations can be studied variationally. There are essentially two types of excitations: (i) solitons connecting the two degenerate spontaneously dimerized ground states (exactly speaking, soliton pairs), and (ii) magnons which are localized excitations and can be viewed as soliton-antisoliton bound states. The variational results give an upper bound for the dispersion at the dimer line where the ground state is exact, and one may expect that they provide a reasonably good approximation in the vicinity of this line.

### A Solitons

For a single soliton excitations one can use the variational ansatz of the type introduced in Ref. :

$$|\alpha ;p^{\sigma \tau }=\underset{n}{}\left\{e^{ipx_n^L}|L_n^{\sigma \tau }+\alpha e^{ipx_n^R}|R_n^{\sigma \tau }\right\}.$$
(4)
Here the states $`|L_n^{\sigma \tau }`$, $`|R_n^{\sigma \tau }`$ are shown schematically in Fig. 1, and $`x_n^L=2n+\frac{1}{2}`$, $`x_n^R=2n+\frac{3}{2}`$ are the effective “center of mass” coordinates of $`|L`$ and $`|R`$ states, respectively, and $`\alpha `$ is a variational parameter. This soliton state has the quantum numbers $`S^z=\sigma `$, $`\tau ^z=\tau `$ and may be combined in pair with the corresponding antisoliton state to get states with $`(S,\tau ^z)=(1,\pm 1)`$, $`(0,0)`$, $`(1,0)`$ or $`(0,\pm 1)`$. Since there are no diagonal exchange interactions in our “ladder”, the variational energy of a single soliton does not depend on $`(\sigma ,\tau )`$ so that we omit the corresponding index in what follows.

It is easy to see that there exist a couple of useful symmetry properties. Upon the operation of interchanging upper and lower leg $`(S\tau )`$ the $`R`$ and $`L`$ soliton states transform as $`|L_n|R_{n1/2}`$, $`|R_n|L_{n+1/2}`$, so that

$$|\alpha ;p\alpha |(\alpha )^1;p.$$
(5)
That means that if the variational energy of the state (4) calculated with the Hamiltonian (3) has a minimum at a certain $`\alpha =\alpha _0`$, then the energy of the same state calculated with the Hamiltonian with interchanged $`S`$ and $`\tau `$ operators should have a minimum at $`\alpha ^{}=1/\alpha `$. An important consequence is that in the symmetric case $`\epsilon =0`$, $`J_s=J_\tau `$ the minimum should always occur at $`\alpha =\pm 1`$. (It should be remarked that this property was missed in our earlier treatment; the variational parameter $`\zeta \alpha e^{ip}`$ used in Ref.  was erroneously assumed to be real, which forced the minimum to be at $`\zeta =\pm 1`$).

Upon a reflection $`xx`$ the soliton states transform as $`|L_n|R_{Nn}`$, $`|R_n|L_{Nn}`$, so that one has

$`|\alpha ;p`$ $``$ $`\alpha e^{ip(2N+1)}|(\alpha )^1;p`$ (6)
$`=`$ $`\alpha e^{ip(2N+1)}|(\alpha ^{})^1;p^{}.`$ (7)
where the asterisk denotes complex conjugation. The energies of two states with mutually conjugate wave functions should be equal, therefore replacing $`\alpha `$ by $`1/\alpha ^{}`$ should not change the variational energy, which means that the minimum is reached on a subspace with $`|\alpha |^2=1`$. Thus, we can set $`\alpha =e^{i\phi }`$ in our variational ansatz (4).

Essentially simple but tedious calculation of the variational dispersion yields

$`E_{\mathrm{sol}}(p)={\displaystyle \frac{F_0(p)+F_c(p)\mathrm{cos}\phi +F_s(p)\mathrm{sin}\phi }{24(2516\mathrm{cos}^2p)(5+4\mathrm{cos}\phi \mathrm{cos}p)}},`$ (8)
where we have used the notation

$`F_0(p)`$ $`=`$ $`[J_\tau (\epsilon +3)+3J_s](383392\mathrm{cos}(2p))`$ (9)
$`+`$ $`(\epsilon +3)(48\mathrm{cos}(4p)+180\mathrm{cos}(2p)93),`$ (10)
$`F_c(p)`$ $`=`$ $`4\mathrm{cos}p(3332\mathrm{cos}2p)\{(2J_\tau 3)(\epsilon +3)+6J_s\},`$ (11)
$`F_s(p)`$ $`=`$ $`20\mathrm{sin}p(2516\mathrm{cos}^2p)\{J_\tau (\epsilon +3)3J_s\}.`$ (12)
It is now easy to minimize the expression (8) in $`\phi `$ and to obtain the final variational soliton energy. The equation $`dE_s/d\phi =0`$ always has two roots, $`\phi _0`$ and $`\pi \phi _0`$, and one of them yields the minimum for $`0<p<\frac{\pi }{2}`$, and the other does so for $`\frac{\pi }{2}<p<\pi `$. In this way, though one cannot see that from Eq. (8) directly, the final minimized soliton energy $`E_{\mathrm{sol}}^{\mathrm{min}}(p)`$ is symmetric with respect to the shift $`pp+\pi `$, in full agreement with the fact that the Brillouin zone in the dimerized phase is halved.

For certain special cases the variational soliton dispersion can be sufficiently simplified. It is easy to see that under the following condition of “generalized symmetry”:

$$J_s=J(1+\epsilon /3),J_\tau =J$$
(13)
the function $`F_s(p)`$ identically vanishes, so that the minimum of (8) is reached at $`\phi =0`$ or $`\pi `$ (i.e., $`\alpha =\pm 1`$), and one obtains

$`E_{\mathrm{sol}}^{\mathrm{sym}}`$ $`=`$ $`(3+\epsilon )/24`$ (14)
$`\times `$ $`{\displaystyle \frac{62J+12\mathrm{cos}2p+3(64J+12)|\mathrm{cos}p|}{54|\mathrm{cos}p|}}.`$ (15)
At the dimer line, i.e., for $`J=\frac{3}{4}`$, the above expression simplifies further to

$$E_{\mathrm{sol}}^{\mathrm{dimer}}=\frac{3+\epsilon }{16}(54|\mathrm{cos}p|).$$
(16)
The energy gap at the dimer line occurs always at $`p=0`$, $`\pi `$. One can check that the gap which follows from (16) coincides with the result given in Ref. ; however, for general $`p`$ the formula (16) yields the lower energy and is in addition much simpler than the corresponding expression in Ref. .

The energy of the continuum of soliton-antisoliton scattering states is given by

$`E_{pair}(k,q)=E_s(k/2+q)+E_a(k/2q),`$

where $`k`$ and $`q`$ are respectively the total and relative momentum. All momenta here are defined for the original non-dimerized chain, and the Brillouin zone is in fact halved, so that $`k[\pi /2,\pi /2]`$.

### B Magnons

We understand “magnons” as localized excitations which do not break the dimer order (in contrast to solitons). In a crude approximation, one could imagine that a magnon state corresponds to exciting one of the singlet bonds on the $`S`$ or $`\tau `$ chains to a triplet and making it propagate along the chain. One may speak about $`S`$\- or $`\tau `$-magnons, respectively. In turns out, however, that such a crude single-mode approximation overestimates considerably the magnon energy. A simple reason is that a state $`|t_{n\mu }`$ shown in Fig. 2 is not an eigenstate of the Hamiltonian even at the dimer line, and under the action of the Hamiltonian it gets mixed with other states. In the vicinity of the dimer line the largest contribution to the first-order matrix elements (i.e., of the $`\psi |\widehat{H}|t_{n\mu }`$ type) comes from the action of the local Hamiltonian $`\widehat{h}_n`$ coupling the sites $`(2n+1)`$ and $`(2n+2)`$ (see Fig. 2). The other first-order contributions are proportional to small parameters $`\lambda _\tau =J_\tau \frac{3}{4}`$, $`\lambda _s=J_s\frac{3+\epsilon }{4}`$, measuring the deviations from the dimer line, and we neglect them for the sake of consistency (see below). Thus, it looks reasonable to take a linear combination of the two orthogonal states $`|t_{n\mu }`$ and

$$|\mathrm{\Delta }^{(1)}t_{n\mu }=\widehat{h}_n|t_{n\mu }t_{n\mu }|\widehat{h}_n|t_{n\mu }t_{n\mu }|$$
(17)
(see Fig. 2) as a simplest variational ansatz. The magnon state with momentum $`k`$ can be written as

$$|m;k=\underset{n}{}e^{ik(2n+1)}\left\{\mathrm{cos}\chi |t_{n\mu }+\mathrm{sin}\chi |\mathrm{\Delta }^{(1)}t_{n\mu }\right\}.$$
(18)
where $`|\mathrm{\Delta }^{(1)}t_{n\mu }`$ is assumed to be normalized, and $`\chi `$ is a variational parameter. After minimization in $`\chi `$, the magnon dispersions take the following form:

$`E_S(k)`$ $`=`$ $`e_{0,S}+(1/8)(3+\epsilon 4J_s)\mathrm{cos}(2k),`$ (19)
$`E_\tau ^0(k)`$ $`=`$ $`e_{0,\tau }^0+(1/8)(1+\epsilon )(34J_\tau )\mathrm{cos}(2k),`$ (20)
$`E_\tau ^{\pm 1}(k)`$ $`=`$ $`e_{0,\tau }^{\pm 1}+(1/8)(34J_\tau )\mathrm{cos}(2k),`$ (21)
where the superscripts for $`\tau `$-magnons indicate the value of $`\tau ^z`$. One may observe that the inclusion of the above-mentioned contributions to $`|\mathrm{\Delta }^{(1)}t_{n\mu }`$ which are proportional to $`\lambda _s`$, $`\lambda _\tau `$ would lead to corrections of the order $`O(\lambda _{s,\tau }^2)`$ in both the bandwidth and the self-energy of magnons. To be consistent, we have to omit those contributions, since we have already implicitly neglected the corrections to the ground state energy which are also $`O(\lambda _{s,\tau }^2)`$ in the leading order. (This does not mean, of course, that our calculation will have the accuracy $`O(\lambda _{s,\tau }^2)`$, since the contributions from the higher-order matrix elements of the type $`\psi |(\widehat{H})^m|t_{n\mu }`$ do not contain any small parameter. The ansatz (18) should be considered as purely variational, sort of “one-step Lanczos” trick).

The self-energies for $`\tau `$-magnons can be compactly written as

$`e_{0,\tau }^\mu `$ $`=`$ $`3J_s/4+J_\tau `$ (22)
$``$ $`{\displaystyle \frac{(1+\eta \epsilon )+2\left\{(1+\eta \epsilon )^2+48J_s^2\right\}^{1/2}}{16}},`$ (23)
where $`\eta =1`$ for $`\mu =\pm 1`$ and $`\eta =1`$ for $`\mu =0`$. The expression for the self-energy of an $`S`$-magnon is somewhat cumbersome due to the fact that the state $`|\mathrm{\Delta }^{(1)}t_{n\mu }`$ in case of the $`S`$-magnons does not have a certain value of the total orbital momentum $`T`$ and is a mixture of states with $`T=0`$ and $`2`$:

$`e_{0,S}`$ $`=`$ $`J_s+J_\tau {\displaystyle \frac{4J_\tau (3\epsilon )+3(1+\epsilon )+\sqrt{Z}}{16(\epsilon ^2+2\epsilon +3)}},`$ (24)
$`Z`$ $`=`$ $`16J_\tau ^2(\epsilon ^2+2\epsilon +2)[\epsilon ^4+4\epsilon ^3+27(\epsilon ^2+2\epsilon +2)]`$ (25)
$`+`$ $`8J_\tau \epsilon ^2(\epsilon ^4+6\epsilon ^3+21\epsilon ^2+32\epsilon +24)`$ (26)
$`+`$ $`\epsilon ^6+6\epsilon ^5+21\epsilon ^4+44\epsilon ^3+36(2\epsilon ^2+2\epsilon +1)`$ (27)
One can observe that the ansatz for a magnon state which we have used above is in fact equivalent to the simplest possible ansatz for a bound state of soliton and antisoliton, where only the configurations of the type

$`|L_n|\overline{L}_{n+\frac{1}{2}},|R_n|\overline{R}_{n+\frac{1}{2}},|L_n|\overline{R}_{n+\frac{1}{2}},|R_n|\overline{L}_{n+\frac{1}{2}}`$

with soliton and antisoliton being nearest neighbors are taken into account. Here the bars denote antisoliton states; note that for fixed boundary conditions, solitons and antisolitons live on different sublattices, e.g. the free spin in the $`𝝉`$ chain (see Fig. 1) is placed on the odd (even) rung for solitons (antisolitons), respectively, which is reflected in the above notation.

This indicates that magnons should be viewed as bound soliton-antisoliton states, rather than as independent excitations. Our simple ansatz can capture only the short-range physics of the tightly bound states with quantum numbers $`(S,\tau ^z)=(1,0)`$, $`(0,0)`$, and $`(0,\pm 1)`$; it will miss, however, possibly existing higher-lying loosely bound states as well as states with $`(S,\tau ^z)=(1,\pm 1)`$.

It is interesting to remark that solitons themselves can be viewed as bound states of two spinons belonging to the $`S`$ and $`\tau `$ chains, respectively, which should unbind when one moves towards the exactly solvable $`SU(4)`$ point. Therefore on the way from the $`SU(4)`$ point to the dimer line there should exist an interesting hierarchy of bound states, which would be an interesting subject of the future studies.

### C Crossover between solitons and magnons

The minimum of the lower boundary of the soliton continuum is at $`k=0`$ in the immediate vicinity of the dimer line. When $`J_s`$, $`J_\tau `$ decrease below certain values, the second equivalent minimum splits off and starts to move towards $`\pi /2`$. The magnon gap always lies either at $`k=\pi /2`$ or at $`k=0`$. Fig. 3 illustrates the behavior of the gap along the “generalized symmetric line” $`J_s=J(1+\epsilon /3)`$, $`J_\tau =J`$ for a few choices of $`\epsilon `$.

One can observe the following important feature: the solitons are the lowest excitations only in a limited region of the phase space. Even at the dimer line line solitons cease to be the lowest excited states if the anisotropy is sufficiently strong: the $`\tau ^\pm `$-magnons become the lowest excitations for $`\epsilon <59/33`$ and soften at the phase transition point $`\epsilon =2`$ (respectively, $`\tau ^0`$-magnons become the lowest excitations for $`\epsilon >59/7`$ and the corresponding gap saturates at the finite value of $`\frac{3}{4}`$ at $`\epsilon +\mathrm{}`$). One can see that whenever the variational gap closes while moving from the dimer line, the magnon gap touches zero first, so that the phase transitions are always governed by magnons.

Three branches of $`\tau `$-magnons are degenerate at $`\epsilon =0`$, and for $`\epsilon >0`$ ($`\epsilon <0`$) the $`\tau ^0`$-magnons ($`\tau ^\pm `$-magnons) become lower in energy, respectively; the $`S`$-branch lies always in between the $`\tau ^\pm `$ and $`\tau ^0`$ branches. In this connection, we would like to make a remark concerning Fig. 4 of Ref.  which shows the behavior of the $`S`$-magnon gap (the lowest excitation in the $`(S^z,\tau ^z)=(1,0)`$ sector) across one of the phase boundaries in the model (3) with $`\epsilon =1`$. We think the behavior of $`S`$-gap is not relevant as a proof of the first order transition, because another excitation, namely that of the $`(0,1)`$ sector, should be responsible for the transition.

Away from the dimer line our description should become progressively worse since the pure dimer state is too far from the correct ground state. Thus the present approach cannot be used for estimation of the shape of the phase boundaries, except when they become close to the dimer line (see below).

Another observation is that though the gap of the soliton continuum can be at incommensurate $`k`$, the lowest excited state almost always (actually except for very narrow regions, whose existence may well be an artifact of our approximation) occurs at the commensurate values of $`k=0`$ or $`\pi /2`$ (see Fig. 3), which is consistent with the results of Refs. . However, the contribution from solitons can be observed in the dynamic structure factor $`S(k,\omega )`$. Unfortunately, the only existing dynamical DMRG study  provides data only for one point in the dimerized phase (namely, for the point $`J_s=\frac{3}{4}`$, $`J_\tau =\frac{3}{4}`$, $`\epsilon =0`$) where no incommensurability is expected.

One can see that if $`\epsilon `$ tends towards the multicritical points $`\epsilon =2`$ or $`\epsilon =+\mathrm{}`$, then there are two transitions along the generalized symmetry line in the vicinity of the point $`J=\frac{3}{4}`$ belonging to the dimer line: one at $`J<\frac{3}{4}`$, and another at $`J>\frac{3}{4}`$. This indicates that if the anisotropy exceeds a certain critical value, an additional transition boundary should appear in the ($`J_\tau `$, $`J_s`$) plane. Note that in this case the variational phase boundaries lie close to the dimer line, so that this statement can be considered reliable.

### D Comparison with numerical results

We have studied the spectrum of low-lying excitations of the model (3) by means of the exact diagonalization (Lanczos) method, using a chain consisting of 12 sites (i.e. having 12 spin and 12 pseudospin degrees of freedom).

Fig. 4 illustrates the behavior of the elementary excitations spectrum for a few choices of $`J_s`$, $`J_\tau `$, $`\epsilon `$, in comparison with the variational results. For each value of the momentum $`k`$, the energies of four lowest-lying eigenstates in each of the subspaces with $`(S^z,\tau ^z)=(0,0)`$, $`(1,1)`$, $`(1,0)`$, and $`(0,1)`$ were calculated. (The subspaces are denoted in Fig. 4 respectively with circles, diamonds, down and up triangles). One can see that our simple variational ansatz is in a qualitative agreement with the Lanczos data for the lowest excited states in the vicinity of the dimer line. It explains the main feature of the Lanczos data, namely, the rather flat dispersions observed nearly everywhere in the dimerized phase. One can also observe the tendency of solitons to have a spectrum with the minimum at an incommensurate value of $`k`$ when $`J_s`$, $`J_\tau `$ decrease below certain value (see Fig. 4b).

Identification of different types of excitations can be done as follows: Free soliton-antisoliton pairs can have any of the four above-mentioned $`(S^z,\tau ^z)`$ combinations, whose energy would be exactly degenerate in the infinite system; that allows one to identify the soliton states in Fig. 4a-d as highly degenerate (in reality, nearly degenerate) points where several different symbols are sitting on top of each other. $`S`$-magnons appear as degenerate combination of circle and down triangle, $`\tau ^0`$-magnons should be circles, and $`\tau ^+`$-magnons must be up triangles.

The agreement concerning the higher-lying excitations is considerably worse. Only for $`\epsilon =0`$ (see Fig. 4a,b) one can see the isolated magnon branch; for high anisotropies in the vicinity of the multicritical point (Fig. 4e,f) a large number of states appears immediately above the lowest-lying states appears, which may mean that actually many bound states split from the soliton continuum, and we are able to capture only the lowest-lying, most tightly bound states due to the short-range nature of the magnon ansatz (18).

## IV Doped spin-orbital model: charged solitons and their bound states

It is interesting to study the behavior of the spin-orbital model doped with holes (i.e., away from the quarter filling). In this case one has to consider the corresponding $`t`$-$`J`$ model, adding the appropriate hopping term to the Hamiltonian $`H(J_s,J_\tau ,\epsilon )`$ determined by Eq. (3):

$$\widehat{H}=\widehat{H}(J_s,J_\tau ,\epsilon )+t\underset{l\sigma \tau }{}(c_{l,\sigma \tau }^{}c_{l+1,\sigma \tau }^{}+\text{h.c.}).$$
(28)
Here $`c_{l,\sigma \tau }`$ are the electron Fermi-operators at the site $`l`$, $`\sigma `$ and $`\tau `$ are the corresponding spin and orbital quantum numbers, and the condition of no double occupancy is implicitly assumed. We also assume, in a similar fashion as in the course of derivation of (1), that hopping exists only between the nearest neighbors with the same type of orbitals; the hopping amplitude $`t`$ is chosen to be negative for the sake of definiteness. One may also wish to add the nearest-neighbor Coulomb repulsion term

$`\widehat{H}\widehat{H}+V{\displaystyle \underset{l}{}}\widehat{n}_l\widehat{n}_{l+1},`$

where $`\widehat{n}_l=_{\sigma \tau }c_{l,\sigma \tau }^{}c_{l,\sigma \tau }^{}`$.

Consider a single hole in this system. If one is on the dimer line $`J_\tau =3/4`$, $`J_s=(3+\epsilon )/4`$, where the checkerboard-type singlet state is exact, it is easy to show that the hole state $`|h;n`$ shown schematically in Fig. 5 gets hopped to the neighboring site without disturbing the structure of the ground state, so that the corresponding translational invariant state with momentum $`k`$ is an exact eigenstate with the energy

$$E_{1h}(k)=2e_0+2t\mathrm{cos}(k),e_0\frac{3(3+\epsilon )}{16}.$$
(29)
Such a state connects two regions with different dimerization pattern, so that it can be viewed as a “charged soliton”, or a hole bound on a soliton. One can straightforwardly see that a hole state which does not disturb the dimer order is not able to move and has to decay into a “charged soliton” and an antisoliton.

The two-hole problem can be also easily studied. The general two-hole state with a total momentum $`k`$ can be written as

$$\mathrm{\Psi }(k)=\underset{n_1<n_2}{}e^{ik(n_1+n_2)/2}f(n_2n_1)|n_1,n_2,$$
(30)
and from the Schrödinger equation one obtains the following system of equations for the function $`f(r)`$:

$`2t\mathrm{cos}{\displaystyle \frac{k}{2}}[f(r+1)+f(r1)]=(E_{2h}4e_0)f(r),r2,`$ (31)
$`2t\mathrm{cos}{\displaystyle \frac{k}{2}}f(2)=(E_{2h}3e_0V)f(1).`$ (32)
Setting $`f(r)=Ae^{iqr}+Be^{iqr}`$, with the parameter $`q`$ (having the meaning of the relative momentum) being real, one obtains the continuum of scattering states with the energy

$$E_{2h}(k,q)=4e_0+4t\mathrm{cos}(k/2)\mathrm{cos}(q/2).$$
(33)
There is also a bound state solution with $`f(r)=f_0e^{\kappa r}`$; the parameter $`\kappa >0`$ is determined by the equation

$$e^\kappa =\frac{Ve_0}{2t\mathrm{cos}(k/2)}.$$
(34)
This solution describes a localized state with the energy

$$E_{2hb}^b(k)=3e_0+V\frac{2t^2}{e_0V}(1+\mathrm{cos}k).$$
(35)
The spectrum of the 2-hole problem is schematically shown in Fig. 6. One can see that for $`V<e_02|t|`$ the energy of this state lies below the lower boundary of the two-hole continuum in the whole range of $`k`$. When $`e_0<V<e_02|t|`$, the bound state exists only in a certain interval $`[k_0,\pi ]`$ where $`k_0`$ is determined from (34) by demanding that $`\kappa `$ is positive; $`k_0`$ tends to $`\pi `$ as $`V`$ tends to $`e_0`$. For $`V>e_0`$ the energy of the localized state (35) lies above the upper boundary of the continuum, so that the states are “antibound.” Thus, the “charged solitons,” which are themselves bound states of a hole and a soliton, tend, in their turn, to form bound states, if the condition $`e_0>V`$ is satisfied.

One can observe that the general $`N`$-hole problem is equivalent to that of the $`XXZ`$ ferromagnet with $`J_z=e_0V`$ and $`J_{xy}=2t`$, in external magnetic field $`h=e_0+V`$ along the $`z`$ axis, and thus can be treated by means of the Bethe ansatz. In fact, one can use the known result for the energy of $`N`$-spin bound complex  to translate it into the corresponding formula for our problem:

$`E_N(k)`$ $`=`$ $`(e_0+V)N+\left\{(e_0V)^24t^2\right\}^{1/2}`$ (36)
$`\times `$ $`\left\{\mathrm{tanh}(N/N_0)+{\displaystyle \frac{2}{\mathrm{sinh}(2N/N_0)}}\mathrm{sin}^2(k/2)\right\},`$ (37)
where $`N_0`$ is defined by the relation

$`\mathrm{cosh}(2/N_0)=(Ve_0)/2t.`$

One can straightforwardly check that (35) is the particular case of (36) for $`N=2`$. For large $`N`$ such a $`N`$-hole complex corresponds to a classical soliton state in the $`XXZ`$ ferromagnetic chain.  The expression (36) is valid for $`|Ve_0|>2|t|`$. It is known  that for $`|Ve_0|<2|t|`$ only bound (respectively antibound) states with $`N=2`$ exist, and the energies of states with larger $`N`$ fall into the continuum. For $`Ve_0<2|t|`$ the $`N`$-hole bound complex is stable against the decay into smaller composite particles,  i.e.,

$`E_N(k)<\mathrm{min}\{E_{NM}(kk^{})+E_M(k^{})\},`$

and the opposite inequality (meaning respectively instability) is fulfilled if $`Ve_0>2|t|`$. Thus, for $`Ve_0<2|t|`$ the ground state of the model (28) at finite doping represents a condensate in which all holes are bound in a single “particle” with zero total momentum; in the intermediate regime $`2|t|<Ve_0<2|t|`$ only the two-hole bound states are stable, and for $`Ve_0>2|t|`$ all bound states become unstable.

## V Summary

We have studied, by means of the combination of variational and exact results, as well as numerically, the interplay between different types of elementary excitations in the model of an anisotropic spin-orbital chain. We analyze the vicinity of the special “dimer line” in the phase space of the model where the spontaneously dimerized ground state is known exactly. Solitons in the dimer order are shown to be the lowest excitations only in a narrow region of the phase diagram, and the phase transitions are always governed by magnon-type excitations which can be viewed as soliton-antisoliton bound states. It is suggested that an additional phase boundary appears when the anisotropy exceeds certain critical value.

We have considered also a $`t`$-$`J`$ version of the doped spin-orbital model. The elementary charge excitation is shown to be a “charged soliton”, or a soliton bound on a hole. We have studied the conditions under which the “charged solitons” tend to form bound states; the exact expression for the energy of the $`N`$-hole bound complex, based on the known results from the Bethe Ansatz, is presented.

###### Acknowledgements.

This work was initiated during the “Quantum Magnetism” conference at the Institute for Theoretical Physics, Santa Barbara, and HJM wishes to acknowledge the hospitality and support of the ITP during this conference. AK acknowledges the hospitality of Hannover Institute for Theoretical Physics. This work was supported by the German Federal Ministry for Research and Technology (BMBFT) under the contract 03MI5HAN5.


format: text

text:
# Microscopic Theory of Protein Folding Rates.I: Fine Structure of the Free Energy Profile and Folding Routes from a Variational Approach

## I Introduction

Considerable progress has been made in describing protein folding using equilibrium and nonequilibrium statistical mechanics, but a complete formal microscopic kinetic theory has only been sketched. The primary novel features of the modern theory of folding revolve around two somewhat different themes  — the glassy dynamics expected for most heteropolymers whose sequence is chosen at random , and the organized dynamics expected for proteins selected by evolution to fold quickly on a funneled energy landscape . The establishment, through selection, of a funneled landscape entails minimizing frustration  — the conflict between different energy contributions in a random sequence. Once a funneled landscape is established by selection, however, it is the interplay between entropy and guiding energies of the funnel that figure most prominently in determining the observed kinetics. If we neglect sidechain and solvent degrees of freedom, the entropy depends crucially on the polymer physics of the protein chain. The statistical mechanical theory of both the stable states and some transition states has already been outlined . In this paper and its companion  we show how the existing framework can be extended to yield a complete microscopic theory of folding rates of completely minimally frustrated proteins. Microscopic calculations of transition state ensembles, activation free energies, and dynamical pre-factors involving chain motions can be obtained. A brief report on the early progress of this work has already appeared , but here we fill in the details and also explore some interesting polymer physics issues that we only touched on previously.

The folding transition can be considered as a finite-size phase transition involving two or more stable phases: one a high entropy denatured state with little structure, although perhaps collapsed; and the other the low entropy folded state . A quantitative way to distinguish among these phases is through the magnitude of the fluctuations of each monomer about its average position. In the denatured states, the protein explores many conformations, and because there is little well defined structure, the fluctuations of a residue about any particular position are relatively large. In the folded state, the conformations are much more restricted and can be described as relatively small fluctuations about the localized positions of the average native structure . These small amplitude deviations are reflected in part by the “temperature factors” (or Debye-Waller factors) which can be experimentally obtained from fitting X-ray crystallography data to a model structure that allows for these fluctuations.

The qualitative difference between the liquid and solid phases in an ordinary first order crystallization transition can also be described by fluctuation magnitudes. The thermodynamics and to some extent the kinetics of first order phase transitions can be studied using a free energy density functional . To analyze the kinetics of the transition one introduces an approximate density profile that is able to interpolate between the two phases, modeling the formation of crystallites or droplets. The stable phases which are minima of the functional and the critical nucleus which is a saddle-point of the functional can then be described through the variational parameters of the trial density. From the point of view of structure, the folded state of the protein is similar to one specific minimum of an amorphous solid in so far as it is not infinitely periodic. One major difference between most inorganic solids and proteins is that proteins are polymers where topology is important; the residues not only have a chemical identity but also a definite sequence. Accordingly, the uniqueness of the folded structure refers not only to the 3-dimensional shape, but also to the specific residues localized at the coordinates of the native structure. This suggests that in order to apply the free energy functional formalism to study the folding of a particular protein, one should use order parameters that are local in sequence.

An approach to folding kinetics based on a “site-resolved” protein folding free energy functional was presented in Ref. jjp:st:pgw:98 . In this approach, a variational free energy surface is introduced directly through a reference Hamiltonian, which provides a good way of approximating the density. The resulting free energy surface is parameterized by the fluctuations of each residue about the average, folded conformation. When put into practice the scheme is very similar in spirit to the density functional calculations described above. The variational approach has also been used to characterize the protein folding phase diagram  as well as to study folding nucleation  without attending to the specifics of a given native protein structure. In order to explore general issues of the role of the polymer chain statistics in folding, we specialize our calculation to study a particular protein that folds to a known structure. The example we choose, the $`\lambda `$-repressor protein, has been much studied in the laboratory  but we will explore how its folding routes and free energy profile are changed upon varying polymer statistics sometimes using “unphysical” values of the backbone parameters in order to gain insight.

## II Gaussian Model for a Stiff Chain

Consider the Gaussian approximation to the probability density for the $`n`$ monomer positions $`\{𝐫_i\}`$ of a polymer chain

$$\mathrm{\Psi }[\{𝐫\}]\mathrm{exp}\left[\frac{3}{2a^2}\underset{ij}{}𝐫_i\mathrm{\Gamma }_{ij}𝐫_j\right],$$
(1)
where $`a`$ is a microscopic length scale taken to be the mean square distance between adjacent monomers, and we have assumed the mean position vanishes. The correlations of monomer positions in Eq.(1) are given by $`𝐫_i𝐫_j/a^2=[\mathrm{\Gamma }^1]_{ij}`$.

Different choices of the correlations result in different Gaussian models for the polymer backbone. It is most natural to define the model chain in terms of the correlations between the $`n1`$ bond vectors: $`𝐚_i=𝐫_{i+1}𝐫_i`$. Denoting correlations between bond vectors by $`𝐚_i𝐚_j/a^2=[\mathrm{\Gamma }^{(\mathrm{b})1}]_{ij}`$, the bond correlations and the monomer position correlations are related by

$$\mathrm{\Gamma }=M^T\mathrm{\Gamma }^{(\mathrm{b})}M,$$
(2)
where $`M`$ is the $`(n1)\times (n)`$ nearest neighbor difference matrix

$$M=\left[\begin{array}{ccccc}1& 1& & \mathrm{}& 0\\ & 1& 1& & \mathrm{}\\ \mathrm{}& & \mathrm{}& \mathrm{}& \\ 0& \mathrm{}& & 1& 1\end{array}\right]$$
(3)
A simple way to account for chain stiffness is to assume there is a fixed angle $`\theta `$ between adjacent bonds. This freely rotating chain model  has monomer correlations that decay exponentially as $`𝐚_i𝐚_{i+l}/a^2=g^l`$, where $`g=\mathrm{cos}\theta `$. Inverting the matrix of bond correlations and using Eq(2) gives the inverse of the monomer correlations

$$\mathrm{\Gamma }=\frac{1g}{1+g}K^R+\frac{g}{1g^2}[K^R]^2\frac{g^2}{1g^2}\mathrm{\Delta },$$
(4)
where $`K^R`$ is the Rouse matrix for a nearest neighbor harmonic chain

$$K^R=\left[\begin{array}{ccccc}1& 1& & \mathrm{}& 0\\ 1& 2& 1& & \mathrm{}\\ & \mathrm{}& \mathrm{}& \mathrm{}& \\ \mathrm{}& & 1& 2& 1\\ 0& \mathrm{}& & 1& 1\end{array}\right],$$
(5)
and $`\mathrm{\Delta }`$ is accounts for the “boundaries” at the end of the chain

$$\mathrm{\Delta }=\left[\begin{array}{cccccc}1& 1& & & \mathrm{}& 0\\ 1& 1& & & & \mathrm{}\\ & & & & & \\ & & & & & \\ & & & & & \\ \mathrm{}& & & & 1& 1\\ 0& \mathrm{}& & & 1& 1\end{array}\right].$$
(6)
In the continuum limit, this Gaussian polymer model gives the familiar form of wormlike chain that restricts the mean separation between adjacent monomers as well as the local curvature of the chain (e.g., see Ref. dt:byh:98  and references therein). In this paper, we use the discrete representation and identify the monomers to be the $`\alpha `$ carbons composing the polypeptide backbone. With this choice, the root mean square separation distance $`a`$ is the typical distance between adjacent $`\alpha `$ carbons $`a3.8`$Å.

The other parameter in the chain model defines the chain stiffness: as $`g0`$, $`\mathrm{\Gamma }=K^R`$ gives the familiar correlations of a flexible chain ($`𝐚_i𝐚_j=a^2\delta _{ij}`$), and as $`g1`$, the correlations correspond to that of a rigid rod ($`𝐚_i𝐚_j=a^2`$). Another measure of the chain stiffness is given by the persistence length $`l_i=_j𝐚_i𝐚_j/a`$ . It is possible to introduce non-uniform stiffness parameters (and hence local persistence lengths) in order to model the different flexibilities of the monomers composing a heteropolymer . For example, the bond correlations can be extracted from the rotational isomeric states model of Flory , or from equilibrium simulations of a detailed polymeric potential . Still more complex models can be used to model explicit chiral helical tendency  through anisotropic Gaussian correlations. For simplicity, we assume here that the chain stiffness is uniform in this paper, so that the persistence length is related to the chain stiffness by $`la/(1g)`$. For proteins, a reasonable value for the chain stiffness is $`g=0.8`$ which corresponds to the persistence length of polyalanine, $`l=5a20`$Å. In this paper we will also compare chains of other uniform stiffnesses as well.

In addition to the Gaussian correlations given by Eq.(1) and Eq.(4) for a stiff chain, we also include a confining potential that controls the overall size of the polymer chain, which is to say it determines the proximity to the chain collapse transition. To model a collapsed stiff chain, we consider the chain Hamiltonian

$$\beta H_{\mathrm{chain}}=\frac{3}{2a^2}\underset{ij}{}𝐫_i\mathrm{\Gamma }_{ij}𝐫_j+\frac{3}{2a^2}B\underset{i}{}𝐫_i^2.$$
(7)
The second term controls the degree of collapse of the chain through a confining potential where the parameter $`B`$ is conjugate to the radius of gyration of the chain. To establish notation for future use, we rewrite Eq.(7) as

$$\beta H_{\mathrm{chain}}=\frac{3}{2a^2}\underset{ij}{}𝐫_i[\mathrm{\Gamma }^{(\mathrm{ch})}]_{ij}𝐫_j,$$
(8)
where

$$\mathrm{\Gamma }_{ij}^{(\mathrm{ch})}=\mathrm{\Gamma }_{ij}+B\delta _{ij}$$
(9)
are the inverse monomer correlations of the collapsed stiff chain.

## III Variational Free Energy Surface

The Hamiltonian for our protein model is

$$H=H_{\mathrm{chain}}+H_{\mathrm{int}},$$
(10)
where $`H_{\mathrm{chain}}`$ is the backbone potential defining the polymeric correlations given in Eq.(8) and $`H_{\mathrm{int}}`$ is the (2-body) interaction potential between distant monomers. The interaction between distant monomers are modeled by a pair potential $`u(r)`$

$$H_{\mathrm{int}}=\underset{[ij]}{}ϵ_{ij}u(|𝐫_i𝐫_j|),$$
(11)
where $`ϵ_{ij}`$, the strength of the interaction, depends on the identity of the residues $`i`$ and $`j`$. The spatial dependence of the interactions between distant monomers consists of an attractive well and a repulsive core. For computational convenience, we approximate the interaction potential as the sum of three Gaussians:

$$u(r)=\underset{k=(\mathrm{s},\mathrm{i},\mathrm{l})}{}\gamma _j\mathrm{exp}\left[\frac{3}{2a^2}\alpha _kr^2\right].$$
(12)
where $`(\alpha _\mathrm{s}>\alpha _\mathrm{i}>\alpha _\mathrm{l})`$ are the ranges of the short-, intermediate-, and long-range interactions, respectively. The intermediate-range term is repulsive $`(\gamma _\mathrm{i}>0)`$ and the long-range term is attractive $`(\gamma _\mathrm{l}<0)`$; the intermediate- and long-ranged potential parameters are chosen so that the sum of these two terms gives a potential well at an appropriate distance for contacts in the native structure. The short-range term is repulsive $`(\gamma _\mathrm{s}>0)`$ and represents the hard core repulsion between residues. (For an example of $`u(r)`$ see Fig. 1)

We approximate the free energy surface of the protein using a reference Hamiltonian that corresponds to a polymer in a non-uniform external field that constrains the monomers to lie near their locations in the native structure $`\{𝐫_i^\mathrm{N}\}`$

$$\beta H_0=\beta H_{\mathrm{chain}}+\frac{3}{2a^2}C_i(𝐫_i𝐫_i^\mathrm{N})^2.$$
(13)
The strengths of the harmonic constraints, $`\{C_i\}`$, are conjugate to the fluctuations of the polymer about each of the native positions. The external constraints in $`H_0`$ influence both the correlations $`G_{ij}`$ and average positions $`\{𝐬_i\}`$ of the monomers composing the reference chain:

$$G_{ij}\delta 𝐫_i\delta 𝐫_j_0/a^2=[\mathrm{\Gamma }^{(0)}]_{ij}^1$$
(14)
$$𝐬_i𝐫_i_0=\underset{j}{}G_{ij}C_j𝐫_j^N,$$
(15)
where $`\delta 𝐫_i`$ is the position of the $`i^{\mathrm{th}}`$ monomer relative to the average

$$\delta 𝐫_i=𝐫_i𝐫_i_0=𝐫_i𝐬_i,$$
(16)
and $`\mathrm{\Gamma }^{(0)}`$ is the matrix of coefficients of the quadratic terms of $`H_0`$

$$\mathrm{\Gamma }_{ij}^{(0)}=\mathrm{\Gamma }_{ij}^{(\mathrm{ch})}+C_i\delta _{ij}.$$
(17)
From the magnitude of these fluctuations, this reference Hamiltonian can distinguish the two stable phases of the protein (as described above): the globule corresponds to large fluctuations (weak constraints) and folded states correspond small fluctuations (strong constraints).

We consider the variational free energy surface parameterized by the constraint parameters $`\{C_i\}`$

$$F[\{C\}]=k_BT\mathrm{log}Z_0+HH_0_0,$$
(18)
where $`Z_0=\mathrm{Tr}\left[e^{\beta H_0}\right]`$ is the partition function of the reference Hamiltonian, and $`\mathrm{}_0=\mathrm{Tr}\left[\mathrm{}e^{\beta H_0}\right]/Z_0`$ denotes the average taken with respect to $`H_0`$. Substituting the expressions for $`H`$ and $`H_0`$ gives the variational free energy $`F=EST`$, where

$$E[\{C\}]=\underset{[ij]}{}ϵ_{ij}u(|𝐫_i𝐫_j|)_0$$
(19)
and

$$S[\{C\}]/k_B=\mathrm{log}Z_0+\frac{3}{2a^2}C_i(𝐫_i𝐫_i^\mathrm{N})^2_0$$
(20)
are the expressions for the energy $`E[\{C\}]`$ and entropy $`S[\{C\}]`$ as functions of the variational constraints.

The averages over $`H_0`$ in Eq.(19) and Eq.(20) can be expressed in terms of $`G`$ and $`\{𝐬_i\}`$, because $`e^{\beta H_0}`$ is a Gaussian distribution. One instructive way to calculate the averages is to introduce approximations to the density of monomer $`i`$, $`\rho _i^1(𝐫)=\delta (𝐫𝐫_i)_0`$,

$$\rho _i^1(𝐫)=\left[\frac{3}{2\pi a^2G_{ii}}\right]^{3/2}\mathrm{exp}\left[\frac{3(𝐫𝐬_i)^2}{2a^2G_{ii}}\right],$$
(21)
and the pair density between $`i`$ and $`j`$, $`\rho _{ij}^2(𝐫)=\delta \left(𝐫(𝐫_i𝐫_j)\right)_0`$,

$$\rho _{ij}^2(𝐫)=\left[\frac{3}{2\pi a^2\delta G_{ij}}\right]^{3/2}\mathrm{exp}\left[\frac{3(𝐫(𝐬_i𝐬_j))^2}{2a^2\delta G_{ij}}\right],$$
(22)
where $`\delta G_{ij}=(\delta 𝐫_i\delta 𝐫_j)^2_0/a^2=G_{ii}+G_{jj}2G_{ij}`$. These densities depend on the constraint parameters $`\{C_i\}`$ through $`G`$ and $`\{𝐬_i\}`$. Averages over $`H_0`$ can be calculated through $`\rho _i^1(𝐫)`$ and $`\rho _{ij}^2(𝐫)`$, for example,

$$u(|𝐫_i𝐫_j|)_0=𝑑𝐫\rho _{ij}^2(𝐫)u(r).$$
(23)
In this way, the variational free energy can be viewed as a density functional with a particular approximation to the density that simultaneously incorporates the polymeric correlations and the monomeric fluctuations about the average positions.

It is straightforward to calculate the entropy and energy in terms of the monomer correlations and mean positions. After some manipulation (see Appendix), the entropy can be written as

$$S[\{C\}]=\frac{3}{2}\mathrm{log}detG\frac{3}{2a^2}\underset{ij}{}s_i\mathrm{\Gamma }_{ij}^{(\mathrm{ch})}s_j+\frac{3}{2}\underset{i}{}C_iG_{ii}.$$
(24)
This expression can be interpreted as follows: the first term is the entropy of the chain due to polymeric fluctuations, the second term is the entropy loss of fixing each monomer to the average positions, and the last term is the entropy of the vibrations about the mean position $`(=(3/2a^2)C_i\delta 𝐫_i^2_0)`$. Similarly, the pair potential can be averaged over $`H_0`$ to give the energy

$$E[\{C\}]=\underset{[ij]}{}ϵ_{ij}u_{ij},$$
(25)
where

$`u_{ij}`$ $`=`$ $`u(|𝐫_{ij}|)_0`$ (26)
$`=`$ $`{\displaystyle \underset{k=(\mathrm{s},\mathrm{i},\mathrm{l})}{}}{\displaystyle \frac{\gamma _k}{(1+\alpha _k\delta G_{ij})^{3/2}}}\mathrm{exp}\left[{\displaystyle \frac{3}{2a^2}}{\displaystyle \frac{\alpha _k(𝐬_i𝐬_𝐣)^2}{1+\alpha _k\delta G_{ij}}}\right].`$ (27)
Finally, we choose to measure the free energy relative to the unconstrained chain

$$\mathrm{\Delta }F[\{C\}]=\mathrm{\Delta }E[\{C\}]T\mathrm{\Delta }S[\{C\}],$$
(28)
where, for example, $`\mathrm{\Delta }F[\{C\}]=F[\{C\}]F[\{C=0\}]`$.

The reference Hamiltonian plays such a prominent role in the variational theory (and in subsequent calculations of folding dynamics) that it warrants further comment. Not surprisingly, other Gaussian models have been previously introduced to model the polymers with fixed contacts or crosslinks. A Hamiltonian of the form

$$\beta H_{\mathrm{h}.\mathrm{c}.}=\frac{3}{2a^2}(𝐫_i𝐫_{i+1})^2+C\underset{[ij]}{}(𝐫_i𝐫_j)^2$$
(29)
has been used to study the thermodynamics  and dynamics  of polymers crosslinked at the sites specified by the pairs $`[ij]`$. This harmonic contact Hamiltonian, $`H_{\mathrm{h}.\mathrm{c}.}`$, has also been used to model the vibrations of folded proteins , with the set $`[ij]`$ limited to the contacts of the native structure, and it was found that the relative magnitudes of monomer fluctuations agree well with measured temperature factors . In contrast to our reference Hamiltonian, $`H_{\mathrm{h}.\mathrm{c}.}`$ is translationally invariant, independent of an explicit native structure. While enforcing this symmetry has some advantages, the potential well is centered at the origin. We note that recently, $`H_{\mathrm{h}.\mathrm{c}.}`$ was used as a mean spherical model for protein folding where a nonlinear constraint was added so that the average monomer positions would lie on the surface of a sphere, preventing the polymer conformation from collapsing to the origin; the minima of $`H_{\mathrm{h}.\mathrm{c}.}`$ with the added imposed condition can yield a meaningful average structure . Nevertheless, the mean locations of the residues are quite distorted. Consequently, $`H_{\mathrm{h}.\mathrm{c}.}`$ is not as well suited as the Hamiltonian we use to describe the protein folding transition where the disordered globule and the structured folded state are separated by a barrier composed of an ensemble of partially ordered conformations. We have, however, investigated $`H_{\mathrm{h}.\mathrm{c}.}`$ as a reference Hamiltonian in the variational context. It gives similar results to those described here.

## IV Order Parameters and Folding Routes

Setting values for the constraint parameters $`\{C\}`$ corresponds to selecting an ensemble of conformations specified by $`\{𝐬\}`$ and $`G`$. The energy of a pair is most stabilizing when the pair density is contained in the potential well, i.e., the mean separation between monomers is within the well and the fluctuations are relatively small. Accompanying this stabilization, however, is the entropy loss of localizing the positions of the pair. In general, when there are many non-zero constraints, the entropy loss due to localization is given by Eq.(24) through the correlations. In other free energy functionals  this entropy loss is estimated but there are difficulties in considering only the entropy loss of individual pairs forming loops because the total entropy loss is inherently nonadditive . The relation of the two approaches is much like the difference between the Thomas-Fermi and Hohenberg-Kohn estimates of the kinetic energy in quantum mechanical density functional theory . The values of the constraint parameters corresponding to the local minima of $`F[\{C_i\}]`$ are a compromise between the energy and entropy decrease of forming contacts. Similarly, the constraints that correspond to saddle-point of $`F[\{C_i\}]`$ also reflect this competition of energy and entropy, because the free energy is minimized in all directions except the unstable mode along which there is a maximum.

The configurational ensemble corresponding to a given set of constraint parameters $`\{C_i\}`$ can also be described by density like order parameters that depend on the local mean square fluctuations. For any set of functions of the chain positions $`\{A_i[\{𝐫\}]\}`$, we can define the order parameter $`A_i[\{C\}]=A_i_0`$, as a function of the constraints. This relationship can be inverted locally provided the Jacobian is nonsingular, $`detJ0`$ with $`J_{ij}=A_i/C_j`$. Since these order parameters are a function of the constraints we can then parameterize the free energy by $`\stackrel{~}{F}[\{𝒜\}]=F[\{C\}]`$ with $`C_i=C_i[\{A\}]`$. For example, the form of $`H_0`$ suggests that the local mean square fluctuations (related to Debye-Waller factors)

$$_i=\delta 𝐫_i^2_0=G_{ii}a^2$$
(30)
are a natural set order parameters for the reference Hamiltonian. (Indeed, this is what motivated our choice of $`H_0`$). In studying the dynamics of barrier crossing in the companion paper, it will prove useful to consider a related but different measure of native similarity. With any identification of the order parameters, we can study the properties of the free energy in $`\{C_i\}`$-space to describe the folding, and then characterize the corresponding ensembles through structural order parameters as equilibrium averages with Hamiltonian $`H_0`$.

We calculate the transition states involved in the folding by searching for saddle-points in $`F[\{C_i\}]`$ using an eigenvector-following algorithm . This algorithm is similar to Newton’s method for optimization, but involves diagonalizing the Hessian matrix, $`^2F/C_iC_j`$, at each iteration. In this routine, the point is updated by stepping in a direction to maximize along the eigenvector with the lowest eigenvalue and minimize along along all others. To find a minimum, a step is taken to minimize along all eigenvectors of the Hessian. In order to use this algorithm, we need to be able to differentiate the free energy with respect to $`\{C_i\}`$, $`_\alpha FF/C_\alpha `$. These derivatives can be easily computed by the chain rule using the elementary derivatives $`_\alpha G_{ij}=G_{i\alpha }G_{\alpha j}`$ and $`_\alpha (\mathrm{log}detG)=G_{\alpha \alpha }`$. The explicit expressions for the derivatives of the energy and entropy are not given here; they are straight-forward to derive, and not very illuminating.

The saddle-points and local minima characterize the average folding routes from this theory. These average pathways are found as follows. The globule and native states are identified by the local minima with the largest and smallest entropy, respectively. These are easy to identify, because the globule is the only stable minimum at high temperature and the native is the only one at low temperature; these minima can be used as the initial guesses for the optimization algorithm for incremental temperature changes until we have these minima at the same temperature. Using linear combinations of these two sets of constraints as initial guesses, we search for a saddle-point. From this saddle-point, we perturb the set of constraints $`\{C_i\}`$ along the unstable eigenvector and use the eigenvector following algorithm with a small step size to find the closest minimum. This gives two local minima, one for each direction on the unstable eigenvector, connected by the saddle-point. This process is repeated until the globule and native state are connected by a series of local minima and saddle-points. We identify this connected sequence as the average folding route, characterizing the transition states and local minima that are important in the folding kinetics. We note that in the example below for the $`\lambda `$-repressor protein, only one folding route was found, but this is not a general result of the theory. The same procedure applied to the SH3 domain has shown that there may be multiple routes (unpublished).

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

In this section, we illustrate the variational theory by studying the folding of a variant of the $`\lambda `$-repressor protein. $`\lambda _{685}`$, is a small (80 residue) protein consisting of five helices in the native structure . This protein folds extremely rapidly following two-state kinetics . From NMR measurements of the folding rate of various mutants, Oas and coworkers concluded that the structure of the transition state consisted mainly of residues in helices H1 and H4 . In Ref. jjp:st:pgw:98  we compared folding routes calculated from the variational theory with these measurements, using a reasonable choice for the persistence length of the polypeptide chain. We investigate here how these results depend on different values of the persistence length (or chain stiffness). These studies allow us to see how some recent simplified approaches to free energy profiles based on assuming complete contiguous sequence folding  become more exact as chain stiffness increases.

### A Model Parameters

To apply the theory, we need to specify the parameters that describe the interaction potential between residues and the polymer chain characteristics. The parameters of the present paper are the same as those chosen in Ref. jjp:st:pgw:98 , but we describe our rationale for this choice in greater detail. For this fast folding protein, we consider a Gō model for the interactions. This means that the sum over residues $`[ij]`$ in Eq.(11) is limited to the set of contacts found in the native structure. This set is defined to be pairs of residues $`(i+4j)`$ that have $`\beta `$ carbons ($`\alpha `$ carbons for glycine) distances within a $`6.5`$Å cutoff in the folded structure. A cutoff between $`68`$Å is commonly used to define contacts in a Gō model, though the precise value is not generally important. We also include in this set residue pairs that are likely to have hydrogen bonds (as determined by the DSSP algorithm ) but fall outside the cutoff. The strength of the interactions for this set depends on the residue identities of the pair. We take the well depths $`ϵ_{ij}`$ to be the magnitude of the Miyazawa-Jernigan energy parameters reported in Ref. sm:rlj:96  in units of $`ϵ_0=k_BT`$.

The parameters for the interaction potential $`u(r)`$ in Eq.(12) are constrained by the $`C_\alpha C_\alpha `$ distances of the set of native contacts. The intermediate- and long-range parameters are chosen so that the sum of these Gaussians has an attractive well that contains all the native contacts and has a minimum value $`u(r^{})=1`$ at the most probable contact distance $`r^{}=1.6a`$. The contact distribution and potential with $`(\gamma _\mathrm{i},\alpha _\mathrm{i};\gamma _\mathrm{l},\alpha _\mathrm{l})`$ = (9,0.54;-6,0.27) are shown in Fig. 1.

The short-range interactions represent the hard core repulsion between the monomers which controls the density of the collapsed polymer. Due to the Gaussian chain approximation, the pair density given by Eq.(22) has non-zero density at short distances, and hence even a de-localized pair has energy contributions from both the repulsive and attractive components of the potential. In the model, the short-range Gaussian amounts to an effective potential that balances the attractive potential for relatively unconstrained polymers. Choosing the repulsion by this criterion is analogous to finding $`\mathrm{\Theta }`$ solvent conditions for the unconstrained polymer (such as the globule). To determine a reasonable repulsive potential for this model, we consider a one-dimensional approximation to the variational free energy by setting all the constraints equal $`C_i=C`$. As can be seen in Fig. 2, the energy as a function of $`C`$ is monotonically stabilizing if $`\gamma _s`$ is small, and has a barrier if $`\gamma _s`$ large. The parameters for the repulsive potential are chosen so that the energy is relatively constant for small values of the constraint parameter: $`(\gamma _\mathrm{s},\alpha _\mathrm{s})=`$ (25,3.0). (This particular value of the strength depends on the width of the short-ranged Gaussian which has chosen somewhat arbitrarily).

The remaining parameter to be specified is $`B`$, the strength of the diagonal confinement term. This confinement parameter is a small constant that effects the fluctuations of unconstrained segments of the chain since the constraint parameters are also diagonal terms in the inverse correlation matrix. Fig. 3 shows the radius of gyration,

$$R_G^2=1/n^2\underset{ij}{}(𝐫_i𝐫_j)^2_0,$$
(31)
evaluated at the native coordinates (dashed) and an unconstrained chain, $`C_i=0`$, (solid) as a function of persistence length for $`B=10^3`$. The unconstrained radius of gyration rises rapidly as the persistence length increases and saturates to a value less than twice the native radius of gyration. The radius of gyration, which in the absence of confinement approaches a large value in the stiff chain limit , is seen to be limited by the confinement term. This tension between local stiffness and confinement is responsible for the plateau in $`R_G`$. Fig. 3 also shows the the radius of gyration evaluated at the constraints corresponding to the globule minimum ($``$). For persistence lengths less than $`l10a`$ the globule constraints lead to a smaller radius of gyration than the free chain value, and for larger persistence length chains the radius of gyrations is somewhat larger. Although it is possible to choose values of $`B`$ in order to set $`R_G`$ of the globule for each chain stiffness, we have chosen to illustrate the effects on the polymer conformations by fixing the confinement and independently varying the chain stiffness.

### B Two-dimensional Illustration

Because the results of calculations in the full variational space (one constraint parameter for each monomer) are somewhat complicated to present, we will illustrate the model in lower dimensions to give the reader an intuitive feeling for the multi-dimensional free energy surface and how folding routes are obtained from the model. We consider a two-dimensional approximation in which we group the protein into two segments and assign a variational spring constant to each group. In this example, monomers with index $`150`$ correspond to $`C_1`$ and the rest correspond to $`C_2`$; These two groups correspond roughly to the helices H1–H3 and H4–H5, respectively. This grouping has a loose correspondence to the folding route obtained from the full variational calculations discussed in the next section.

The free energy surface shown in Fig. 4 has two distinct low energy paths determined by the saddle-points connecting the globule (G) and native (N) states. The average folding routes, as defined by the path from the saddle points to the local minima, are determined by the eigenvector-following algorithm. As can be seen in Fig. 4 these routes are very close to the steepest descents path. Along the two routes from G to N, in Path 1 (dotted line) the constraint $`C_1`$ progressively increase followed by the increase of constraint $`C_2`$, whereas along path 2 (solid line) this order is reversed with $`C_2`$ increasing before $`C_1`$.

The free energy of these paths is plotted parametrically versus the energy ($`\mathrm{\Delta }F(C_1,C_2)`$ vs. $`\mathrm{\Delta }E(C_1,C_2)`$) in Fig. 5, giving a free energy profile where the saddle-points appear as local maxima. Path 2 is the relatively favored route since it has the lower barrier to folding. We note that a rough description of the free energy profile can be represented by connecting the stationary points of $`F(C_1,C_2)`$ by straight lines.

The ensemble of conformations composing the average folding route can be characterized by the magnitude of the fluctuations of each residue $`(_i=\delta 𝐫_i^2_0)`$ for any set of constraints along these paths. Fig. 6 shows the monomer fluctuations evaluated at the constraints corresponding to the globule state and the saddle-points of both paths. These fluctuations give a description of folding consistent with the two-dimensional surface parameterized by the constraints: in path 1 residues in helices H1–H3 become structured at $`\mathrm{TS}_1`$ followed by H4–H5 at $`\mathrm{TS}_2`$, and in path 2 the order is reversed. For a given set of the constraints (i.e., a given saddle-point), the fluctuations are seen to smoothly interpolate between the two groups of monomers. The precise shape of the interface depends on the value of the chain stiffness. A stiffer chain would tend to suppress variations as a function of sequence index resulting in an more gradual interface. As will be seen in the full dimensional calculations, decreasing the chain stiffness allows the magnitude of the fluctuations to change more rapidly between successive monomers, as expected.

### C Fine Structure of the Free Energy Profile: Multiple Transition States and Folding Routes

We follow the same analysis outlined in the two-dimensional illustration to describe the folding paths calculated in the full variational space, but the free energy profile is evaluated only at the saddle-points and local minima rather than as a continuous path.

In this study, we focus on how the folding routes depend on the persistence length of the chain. To put various parameters into context, the homopolymers polyglycine, polyalanine and polyproline have persistence lengths $`l=`$ 6Å $`(2a)`$, 20Å$`(5a)`$, and 220Å $`(60a)`$, respectively. . In the freely rotating chain model $`l=a/(1g)`$, so that these persistence lengths correspond to the chain stiffness parameters $`g=`$ 0.5, 0.8, and 0.98, respectively. Modeling the protein backbone with a single uniform chain stiffness parameter, we take the chain stiffness of polyalanine to be a reasonable value for the protein backbone.

We report the free energy profile at the folding transition temperature, $`T_\mathrm{f}`$. This is the temperature at which the folded and globule ensembles have the same equilibrium probability $`(\mathrm{\Delta }F_\mathrm{G}=\mathrm{\Delta }F_\mathrm{N})`$. Since the entropy and energy loss in these states depend on the chain stiffness, $`T_\mathrm{f}`$ is different for different persistence lengths. (For a flexible chain $`l=2a`$, $`k_BT_\mathrm{f}/ϵ_01.2`$, whereas for a very stiff chain $`l=20a`$, $`k_BT_\mathrm{f}/ϵ_02.2`$). Similarly, the unconstrained ensemble (which we define as the zero of the free energy) is dependent on the persistence length. Consequently, to compare the folding profiles of chains with different persistence lengths, it is convenient to plot the free energy profile relative to the globule free energy against a normalized energy coordinate

$$\overline{E}=(\mathrm{\Delta }E\mathrm{\Delta }E^\mathrm{G})/(\mathrm{\Delta }E^\mathrm{N}\mathrm{\Delta }E^\mathrm{G}),$$
(32)
where $`\mathrm{\Delta }E^\mathrm{G}`$ and $`\mathrm{\Delta }E^\mathrm{N}`$ are the globule and native energy changes, respectively. $`\overline{E}`$ is the fractional stabilization energy and equals 0 at the globule state and 1 at the native state.

Fig. 7 shows the free energy profile for chain stiffness parameters ranging from $`g=0.5`$ to $`g=0.95`$ ($`l=2a`$ to $`l=20a`$). The free energy profile versus energy for a flexible chain ($`l=2a`$) is shown as the solid curve in Fig. 7a. The profile exhibits five transition states (saddle-points) separated by local minima. We specify the transition states sequentially from the globule minimum (G) to the native minimum (N), and the minima by the index of the adjacent transition state (e.g., the minimum between $`\mathrm{TS}_1`$ and $`\mathrm{TS}_2`$ is denoted by $`\mathrm{min}_{12}`$). The profile can be described as a rugged fine structure resulting from the different structural ensembles of the local minima superimposed on a single average free energy barrier. The fine structure on the profile is modest in magnitude amounting to a stabilization of a “high energy intermediate” by at most 1$`k_BT_\mathrm{f}`$. This fine structure should not be confused with the ruggedness due to frustration , which we call “transverse” ruggedness, coming from degrees of freedom different from the one plotted; instead, the fine structure is better described as “longitudinal” ruggedness (along the reaction coordinate) . This is a common feature of many models even those which consider only native contacts  and which therefore have perfectly funnel-like surfaces.

Starting with the most flexible chains shown in Fig. 7a, as the stiffness increases the magnitude of the barrier increases and the the minima separating the transition states become relatively still more shallow (the longitudinal ruggedness diminishes). A flexible chain can take advantage of particularly strong contacts while losing a relatively small amount of entropy to localize the pair. This gives rise to the possibility of relatively stable local minima and lower free energy barriers. Stiffer chains must localize larger segments of the chain (on the order of the persistence length) resulting fewer distinct local minima and larger free energy barriers. For the largest chain stiffness in Fig. 7a ($`l9a`$) the local minima have nearly disappeared, leaving a free energy profile with a single transition state ensemble. The profile for even stiffer chains ($`l=10a`$ to $`l=20a`$) is shown in Fig. 7b. In this parameter range, as the chain stiffness increases further the barrier decreases and the single transition state occurs at a larger fractional stabilization energy (i.e., the energy of transition state becomes closer to the native energy). When the persistence length is large enough, large sections of the chain must be constrained resulting in a more folded transition state. Evidently, for very stiff chains the transition state is similar enough to the native minimum that the barrier decreases.

The barrier height is plotted as a function of persistence length in Fig. 8 (solid line). The maximum barrier height occurs for a persistence length $`l10a`$ with an increase of approximately $`70\%`$ relative to the barrier height for the most flexible chain considered. However, the appropriate energy scale for the transition is $`k_BT_\mathrm{f}`$. In units of $`k_BT_\mathrm{f}`$, the barrier height is relatively constant for a wide range of persistence lengths (dashed line), changing only approximately $`10\%`$ for persistence lengths up to $`l12a`$. The robustness of the barrier height at $`T_\mathrm{f}`$ is interesting since the persistence length of real proteins in the laboratory is not precisely known.

Each ensemble of structures along the average folding route can be characterized by the local temperature factors (Eq.30). The fluctuations corresponding to the transition states for a flexible chain ($`l=2a`$) are shown as the dotted curves in Fig. 9a. The fluctuations corresponding to the local minima are within those of adjacent saddle-points; for clarity, only one local minimum, $`\mathrm{min}_{23}`$, is shown (solid line). The folding route can be characterized by considering the structure of the transition states from the globule to the native structure. The first transition state ensemble $`\mathrm{TS}_1`$ is described by the ordering of helices H4–H5, stabilized by the partial localization of a region of helix H1, while residues in helices H2–H3 remain de-localized. In the subsequent transition states, helix H1 becomes progressively more ordered, while helices H2–H3 continue to have large fluctuations and are the last to order. This scenario for the folding of $`\lambda `$-repressor agrees with the interpretation of kinetic data based on $`\varphi `$-value analysis .

The temperature factors for a larger persistence length ($`l=5a`$) are shown in Fig. 9b. The behavior of the fluctuations describing the transition state ensembles is qualitatively similar to the more flexible chain, though the magnitude of the fluctuations of the disordered regions is larger. Some of the detailed features shown in Fig. 9a have been smoothed out since the chain stiffness suppresses variations along the sequence smaller than the persistence length. For example, the very specific localization of a helix H1 residue shown in the $`\mathrm{min}_{23}`$ curve of the more flexible chain (Fig. 9a) has been broadened to a larger region for the stiffer chain. These differences are rather subtle, but the comparison is useful to illustrate the progression to larger chain stiffnesses. The fluctuations for the single transition state characteristic of larger persistence lengths is shown in Fig. 9c. (Note, here the different curves correspond to different persistence lengths rather than transition states along the same folding route.) While the general shape of the curve is maintained, the magnitude of the fluctuations become increasingly more like those of the native state, and in the largest persistence length considered $`l=20a`$ only the end segments of the chain are significantly disordered.

### D Contiguous Sequence Approximations

Several recent and apparently accurate estimates of folding kinetic parameters have assumed that the transition state ensemble can be described by assigning contiguous segments to be either folded or not and allowing the sequence to be passed into such fully folded or unfolded configurations . As Plotkin et al. argued, such a contiguous sequence approximation (CSA) should apply to late transition states where the entropy is low . How does the polymer’s characteristics determine the quality of this approximation?

Following the contiguous sequence approach, to simplify the problem we reduce the number of states in the variational theory’s description of the protein ensembles by restricting the conformations to only those in which the structured residues are fully native and contiguous in sequence (single CSA), or alternatively, we consider two contiguous stretches of native residues (double CSA) . For each fixed value of the number of folded residues, $`N_\mathrm{f}`$, we find the minimum free energy configuration (satisfying the single or double contiguous constraint). In this way we can construct a free energy profile as a function of $`N_\mathrm{f}`$. This construction neglects the connectivity of the path since the minimum free energy configuration with $`N_\mathrm{f}`$ folded residues may not be simply related to that with $`N_\mathrm{f}+1`$. This connectivity is an added complication to the approach and can be treated by the methods presented Refs. ea:db:99 . For the purposes of this illustration, we neglect this aspect of the approximation. The approach outlined here is a great simplification of the more complete variational formalism since it avoids the relatively difficult numerical calculation of finding saddle-points as a function of degree of ordering. Still, the number of configurations needed in the double CSA is quite large, but tractable ($`1.7`$ million for $`\lambda `$-repressor). On the other hand, these approximations are rather restrictive since they neglect partial ordering and provide a less microscopic characterization of the folding route.

One issue in comparing the exact variational and the CSA results is how to measure the partial order characterizing the stationary ensembles in a global way to comparable the number of fully native residues in the CSAs. For a given stationary point, we can compare by estimating $`N_\mathrm{f}`$ through the normalized fluctuations

$$\overline{}_i=(_i_i^\mathrm{G})/(_i^\mathrm{N}_i^\mathrm{G}),$$
(33)
where the superscripts G and N denote the fluctuations evaluated at the globule and native state, respectively. As a rough approximation, we define $`N_\mathrm{f}`$ to be the number of residues with $`\overline{}_i0.95`$.

The free energy profile from the variational theory (dashed), the single CSA (long-dashed), and double CSA (solid) are plotted as a function of $`N_\mathrm{f}`$ for three different persistence lengths in Fig. 10. The barrier heights from the simpler approximations are about twice the barriers obtained from the variational calculation for each persistence length. This is consistent with results of the other contiguous sequence approaches, where the barriers are much larger than obtained from simulations of Gō models with pairwise additive forces . In the present comparison, one could also expect this behavior, since a variational Hamiltonian with fewer degrees of freedom generally gives higher barriers (c.f. the barriers in the two-dimensional illustration Fig. 5 with the full dimensional calculation in Fig. 7). Nevertheless, the simpler approximations provide an intuitive explanation that agrees qualitatively with the magnitudes of the barrier heights. Confining the residues to be either folded or unfolded (in contiguous stretches) is responsible for the larger barrier heights, i.e., partial ordering reduces the barrier. This reduction seen in the more accurate calculation is reminiscent of the way wetting between two stable phases in nucleated phase transformations lowers nucleation barriers by reducing surface tension, and is discussed in the context of protein folding in Ref. st:pgw:97b .

Considering the free energy profiles for the two approximations for a flexible chain ($`l=2a`$) plotted in Fig. 10a, the maximum free energy for the single CSA is approximately $`50\%`$ larger than the barrier from the double CSA and occurs at a greater value of $`N_f`$. Again there is no surprise in finding the barrier from the double contiguous approximation is lower, since the single contiguous configurations are a subset of the double contiguous configurations. It is interesting, however, to compare this difference as the chain stiffness increases. For a larger chain stiffness corresponding to polyalanine ($`l=5a`$) shown in Fig. 10b, the difference in the barrier height from the single and double CSAs decreases, with the single CSA barrier approximately $`30\%`$ larger than the double CSA barrier. Increasing the stiffness further ($`l14a`$), the free energy profiles shown in Fig. 10c are more similar, differing by only $`10\%`$. Since the profiles become still more similar as the chain stiffness increases (even though the conformations in the double CSA are less restricted), this comparison suggests that the single CSA more accurately describes stiff chain conformations than it does flexible chain conformations.

To make this connection more precise, we consider which residues are ordered along the folding route in these approximations. In the present treatment, this information can be easily represented by a plot of the folded regions specified by the monomer index as a function of $`N_f`$. The structured parts of the chain are indicated by the shaded regions in Fig. 11. To characterize the structure at the saddle-points of the variational free energy surface we consider the Gaussian measure to the native structure

$`\rho _i`$ $`=`$ $`\mathrm{exp}\left[{\displaystyle \frac{3}{2a^2}}\alpha ^\mathrm{N}(𝐫_i𝐫_i^\mathrm{N})^2\right]_0`$ (34)
$`=`$ $`(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]`$ (35)
This measure of the monomer density relative to the native position is the order parameter employed in the companion paper to study the dynamics of the barrier crossing. The degree of native structure at the transition state can be characterized by the normalized measure

$$\overline{\rho }_i=(\rho _i\rho _i^\mathrm{G})/(\rho _i^\mathrm{N}\rho _i^\mathrm{G})$$
(36)
where the superscripts G and N denote Eq.(34) evaluated at the constraints corresponding to the globule and native states, respectively.

Consider first the folded residues from the double CSA for a flexible chain ($`l=2a`$) shown in Fig. 11a. Consistent with the full variational results, the shaded region clearly shows the structure forms between within helices H4–5 and helix H1. Superimposed on the regions of folded residues (from the double CSA) are $`\rho _i`$ as function of monomer index for the four main transition states of the folding route of a chain of the same stiffness (from the variational theory). The density plots of $`\rho _i`$ indicate that the relatively unfolded regions are indeed more structured than the globule value $`\rho _i=0`$. In particular, at the interface between the folded and unfolded regions in the double CSA, $`\rho _i`$ obtains intermediate values indicating partial ordering. Nevertheless, the folded residues in the double CSA agrees qualitatively with the structure obtained from the variational saddle-points. In contrast, the folded residues in the single CSA do not agree with the variational theory very well at all, as shown in Fig. 11b.

Figs.11c and 11d show structured regions along the folding routes for a greater chain stiffness $`l=5a`$, the value appropriate for polyalanine. For this chain stiffness, the unfolded region between helix H1 and helices H4–H5 in the double CSA closes at a smaller value of $`N_f`$ compared with the more flexible chain. In this sense the structured residues from the two CSAs are in better qualitative agreement, though the discrepancy between the two is still pronounced. This difference between the two approximations is greatly reduced for a chain with a much larger persistence length ($`l14a`$) shown in Figs.11e and 11f. While there are still unfolded regions between folded regions in the double sequence approximation, they only persist for a very limited range of $`N_f`$. For this persistence length, the variational free energy surface has only one transition state as indicated by the density plot of $`\rho _i`$. The structure determined by the saddle-point agrees qualitatively with the structure indicated by double CSA, but also with the single CSA (because the two are similar).

## VI Conclusion

In this paper, we used a variational approach to calculate the free energy profiles and characterize the transition state ensembles applicable when nonnative contacts can be ignored. Using $`\lambda `$-repressor as an illustrative example, we investigated the role of chain stiffness on the fine structure of the free energy profile. We found that increasing the persistence length of the chain tends to smooth the free energy profile, making longitudinal ruggedness less pronounced. The transition state ensemble with very stiff chains was found to be more folded than the ensemble with more flexible chains. These results can be interpreted in terms of the tension between taking full advantage of strong local contacts while still respecting the bending rigidity of the chain. We also found that while the absolute barrier height has a pronounced maximum as a function of persistence length, the barrier scaled by the folding transition temperature $`k_BT_\mathrm{f}`$ was relatively robust over a wide range of persistence lengths.

This study allowed us to investigate the applicability of simpler contiguous sequence approximations proposed recently. Both the free energy profiles and the folded residues along the folding routes suggest that the single CSA more accurately describes stiffer chains. Since the folded residues in the single CSA is roughly independent of the persistence length, the single CSA sequence approximation corresponds to a chain that has a longer effective persistence length than is appropriate for most natural proteins. The double CSA is able to capture the appropriate folded regions near the free energy barrier. On the other hand, the neglect of partial ordering leads to over-estimates of the absolute barrier height for the $`\lambda `$-repressor protein. Nevertheless, these approximations should be good enough to compute the perturbations of the activated free energy found in protein engineering experiments that probe the transition state ensemble ($`\varphi `$-analysis) .

## Acknowledgments

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

## A

In this appendix, we outline a derivation of the expression for the entropy of the constrained chain given in Eq.(24).

The partition function of the harmonically constrained chain is given by

$`Z_0`$ $`=`$ $`\mathrm{exp}\left[A{\displaystyle \underset{i}{}}C_i𝐫_i^\mathrm{N}\right]`$ (A1)
$`\times `$ $`\mathrm{Tr}\mathrm{exp}\left[A{\displaystyle \underset{ij}{}}𝐫_i\mathrm{\Gamma }^{(0)}𝐫_j+2A{\displaystyle \underset{i}{}}C_i\mathrm{r}_i^\mathrm{N}𝐫_i\right],`$ (A2)
with $`A=3/2a^2`$ and $`\mathrm{Tr}`$ denotes $`\mathrm{\Pi }_i𝑑𝐫_i`$. These integrals are easy to evaluate by completing squares

$`Z_0`$ $`=`$ $`(\pi /A)^{3n/2}(detG)^{3/2}`$ (A3)
$`\times `$ $`\mathrm{exp}\left[A{\displaystyle \underset{i}{}}C_i(𝐫_i^\mathrm{N})^2+A{\displaystyle \underset{ij}{}}C_i𝐫_i^\mathrm{N}G_{ij}C_j𝐫_j^\mathrm{N}\right],`$ (A4)
where $`G=\mathrm{\Gamma }^{(\mathrm{ch})1}`$. We note that Eq.(15) can be used to express the second term in the exponent as

$$A\underset{ij}{}C_i𝐫_i^\mathrm{N}G_{ij}C_j𝐫_j^\mathrm{N}=A\underset{i}{}C_i𝐫_i^\mathrm{N}𝐬_i$$
(A5)
in terms of the average monomer positions, $`\{𝐬_i\}`$. From Eq.(20), combining the partition function with

$$A\underset{i}{}C_i(𝐫_i𝐫_i^\mathrm{N})^2_0=A\underset{i}{}C_i(G_{ii}a^2+𝐬_i^2+2𝐬_i𝐫_i^\mathrm{N}+(𝐫_i^\mathrm{N})^2)$$
(A6)
gives the entropy of the constrained chain (ignoring a constant factor)

$`S[\{C\}]={\displaystyle \frac{3}{2}}\mathrm{log}detG`$ $`+`$ $`{\displaystyle \frac{3}{2}}{\displaystyle \underset{i}{}}C_iG_{ii}`$ (A7)
$`+`$ $`A{\displaystyle \underset{i}{}}C_i𝐬_i^2A{\displaystyle \underset{i}{}}C_i𝐫_i^\mathrm{N}𝐬_i`$ (A8)
This expression is easier to interpret after simplifying the last two terms. Inserting the identity, $`G\mathrm{\Gamma }^{(0)}=\mathrm{𝟏}`$, into the last term and using Eq.(15) gives

$`A{\displaystyle \underset{i}{}}C_i𝐫_i^\mathrm{N}𝐬_i`$ $`=`$ $`A{\displaystyle \underset{ikj}{}}C_i𝐫_iG_{ik}\mathrm{\Gamma }_{kj}^{(0)}𝐬_j`$ (A9)
$`=`$ $`A{\displaystyle \underset{kj}{}}𝐬_k\mathrm{\Gamma }_{kj}^{(0)}𝐬_j.`$ (A10)
In this form, the last two terms of Eq.(A7) can be combined

$$A\underset{kj}{}𝐬_k(\mathrm{\Gamma }_{kj}^{(0)}C_k\delta _{kj})𝐬_j=A\underset{kj}{}𝐬_k\mathrm{\Gamma }_{kj}^{(\mathrm{ch})}𝐬_j,$$
(A11)
where we have used the definition of $`\mathrm{\Gamma }^{(0)}`$ given in Eq.(17). Thus, we have

$$S[\{C\}]=\frac{3}{2}\mathrm{log}detG+\frac{3}{2}C_iG_{ii}\frac{3}{2a^2}\underset{kj}{}𝐬_k\mathrm{\Gamma }_{kj}^{(\mathrm{ch})}𝐬_j,$$
(A12)
which is Eq.(24).