text:
# Bibliography

April 15, 2003

Simple Two-State Protein Folding Kinetics Requires

Near-Levinthal Thermodynamic Cooperativity

Hüseyin KAYA and Hue Sun CHAN

Protein Engineering Network of Centres of Excellence (PENCE),

Department of Biochemistry, and Department of Medical Genetics & Microbiology, Faculty of Medicine, University of Toronto, Toronto, Ontario M5S 1A8, Canada

Running title: Origins of Simple Two-State Protein Folding

Key words: calorimetric cooperativity / single-domain proteins / chevron plot /
chevron rollover / non-Arrhenius kinetics / Gō models / contact order /

Corresponding author.
E-mail address of Hue Sun CHAN: chan@arrhenius.med.toronto.edu
Tel: (416)978-2697; Fax: (416)978-8548
Mailing address: Department of Biochemistry, University of Toronto, Medical Sciences Building – 5th Fl., 1 King’s College Circle, Toronto, Ontario M5S 1A8, Canada.

Abstract

Simple two-state folding kinetics of many small single-domain proteins are characterized by chevron plots with linear folding and unfolding arms consistent with an apparent two-state description of equilibrium thermodynamics. This phenomenon is hereby recognized as a nontrivial heteropolymer property capable of providing fundamental insight into protein energetics. Many current protein chain models, including common lattice and continuum Gō models with explicit native biases, fail to reproduce this generic protein property. Here we show that simple two-state kinetics is obtainable from models with a cooperative interplay between core burial and local conformational propensities or an extra strongly favorable energy for the native structure. These predictions suggest that intramolecular recognition in real two-state proteins is more specific than that envisioned by common Gō-like constructs with pairwise additive energies. The many-body interactions in the present kinetically two-state models lead to high thermodynamic cooperativity as measured by their van’t Hoff to calorimetric enthalpy ratios, implying that the native and denatured conformational populations are well separated in enthalpy by a high free energy barrier. It has been observed experimentally that deviations from Arrhenius behavior are often more severe for folding than for unfolding. This asymmetry may be rationalized by one of the present modeling scenarios if the effective many-body cooperative interactions stablizing the native structure against unfolding is less dependent on temperature than the interactions that drive the folding kinetics.

INTRODUCTION

A logical test of any general conception about the driving forces in protein folding is to ascertain whether polymer models incorporating the given idea can predict generic behavior of real proteins.1,2 In these considerations, self-contained polymer models with explicit chain representations3 are of particular importance. Quite obviously, the relationship between model energetics and conformational distribution can only be addressed in a physically plausible manner when chain connectivity and excluded volume are adequately taken into account. Using this analytical framework, we found that even mundane protein properties such as calorimetric two-state cooperativity4-6 and simple two-state folding/unfolding kinetics7,8 are remarkable feats from a polymer standpoint. Simply put, it is nontrivial to construct heteropolymer models with commonly used model interaction schemes to reproduce such behavior. These include popular 2-, 3-, 20-letter models and traditional Gō models (see below).4-8 Generic protein properties thus present severe constraints on modeling. Hence, insight into real protein energetics can be gained by requiring self-contained polymer models to satisfy such constraints.

Motivated by the proposed consistency principle9 or principle of minimal frustration10 for protein energetics, Gō11 and Gō-like models (see, e.g., refs. 8, 12–16 and references therein) have long been used in protein folding investigations. These models postulate that only intrachain interactions found in the native (ground-state) conformation are favorable, all other possible intrachain interactions are assumed to be either neutral or unfavorable. Recently, this native-centric approach to modeling has often been justified as well by the discovery that folding rates of natural small single-domain proteins are well correlated with the contact order17 of their native structures. How well do common Gō models mimic the generic properties of small single-domain proteins? For chain models configured on two-dimensional square lattices, we found that even with their explicit native biases, the common Gō interaction scheme falls far short of producing the type of calorimetric two-state cooperativity observed for many small proteins.4 Three-dimensional Gō-like lattice5-7 and continnum8 models are more proteinlike in this regard, as many of them may be considered calorimetrically cooperative if certain lattitude is allowed for empirical baseline subtractions.5

However, common Gō-like schemes are apparently not capable of producing simple two-state folding/unfolding kinetics with linear chevron plots.18 This we have recently demonstrated in several examples,7,8 including lattice and continuum (off-lattice) models as well as models with a rudimentary implicit-solvent treatment of desolvation barriers.8,19 Thus, the inability of common Gō-like constructs to predict simple two-state folding/unfolding kinetics is not an artifact restricted only to lattice Gō models. Most likely, it is a fundamental problem arising from the additive nature of common Gō-like interaction schemes. Our results indicate strongly that such model interaction schemes — on-lattice and otherwise — afford insufficient cooperativity to capture real two-state protein energetics. Here we address this basic question by constructing and testing novel native-centric lattice models that go beyond common additive Gō-like schemes, with intrachain potentials that can lead to simple two-state folding/unfolding kinetics. The ultimate goal of this line of inquiry, of which the present lattice exercise is only a first step, is to decipher the many-body cooperative interactions underlying the behavior of small single-domain proteins.

MODELING A COOPERATIVE INTERPLAY BETWEEN LOCAL
CONFORMATIONAL PREFERENCE AND PROTEIN CORE
FORMATION

We first explore several native-centric variants of a 55mer model (Figures 1–4). Their basic features are derived from the original model we put forth recently.6 For the Gō-like constructs studied here, we retain contributions from the term disfavoring the initiation of left-handed helices (equation 1 of ref. 6). Contributions from the 5-letter contact energies are retained for the native contacts in the ground-state conformation, whereas all nonnative contacts are assigned zero energy as in common Gō models. The resulting native-centric model has the energy function

$$E=E_{\mathrm{contact}}^{}+\gamma _{\mathrm{lh}}N_{\mathrm{lh}},$$
$`(1)`$
where the prime superscript on the $`E_{\mathrm{contact}}^{}`$ term indicates that the sum of pairwise 5-letter energies is restricted to native contacts, and the second term on the right disfavors left-handed helices. Here we use the same contact energies and $`\gamma _{\mathrm{lh}}`$ parameter as in ref. 6. The ground-state energy of the present model equals $`36.1`$. We refer to this as model (i).

Next we consider a model that embodies the idea of a cooperative interplay between local conformational preference and the (nonlocal) packing of the protein core.4,6,7 This is motivated by the observation that secondary structure formation in globular proteins is often context dependent, and that short helices are often not stable in isolation but are stable when packed in the core of a protein (see, e.g., ref. 20 and references therein). Here we mimic this effect by assigning a favorable energy $`_{\mathrm{coop}}`$ to each incidence of the conformational situation described in Figure 1, leading to the energy function

$$E=E_{\mathrm{contact}}^{}+\gamma _{\mathrm{lh}}N_{\mathrm{lh}}+_{\mathrm{coop}}h_\mathrm{c},$$
$`(2)`$
where $`h_\mathrm{c}`$ counts the incidences of a fully formed native helix having the cooperative interactions defined by Figure 1, and $`0h_\mathrm{c}4`$ for the present 55mer model. This we refer to as model (ii).

To account for the possibility that optimal packing of the protein core as a whole can impart significant thermodynamic stability to the native structure, we consider yet another model with an extra favorable energy assigned only to the native conformation. The energy function now becomes

$$E=E_{\mathrm{contact}}^{}+\gamma _{\mathrm{lh}}N_{\mathrm{lh}}+_{\mathrm{coop}}h_\mathrm{c}+E_{\mathrm{gs}},$$
$`(3)`$
where the augmented $`E_{\mathrm{gs}}`$ term takes a nonzero favorable value only when the chain is in its unique ground-state conformation. We refer to this as model (iii). We note that the $`_{\mathrm{coop}}`$ and $`E_{\mathrm{gs}}`$ terms introduced in equations 2 and 3 are many-body in nature. Many-body interactions have been investigated in the context of protein folding (see, e.g., refs. 6, 21, 22). However, their relationship with linear chevron plots and simple two-state kinetics has not been much explored.

Figure 2 shows that the many-body cooperative interactions introduced above enhance thermodynamic cooperativity. In calculating the heat capacities of these models, we made the simplifying assumption that the interactions are temperature independent, and set enthalpy equal to the model energy, as in our previous investigations.5-8 For calorimetric two-state behavior, the van’t Hoff to calorimetric enthalpy ratio $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ has to be close to one.4-6 Now the $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratio ($`\kappa _2`$ in ref. 5 without empirical baseline subtraction) equals $`0.804`$ for model (i) which does not contain many-body cooperative interactions. But it is considerably higher at $`0.88`$ and $`0.91`$, respectively, for models (ii) and (iii) with favorable values of $`_{\mathrm{coop}}`$ and $`E_{\mathrm{gs}}`$ (Figure 2, upper panel).*** For every model considered in Figure 2, the $`\kappa _2^{(\mathrm{s})}`$ (ref. 5) value for the $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratio after empricial baseline subtraction equals 1.0. This is not too surprising because the many-body cooperative interactions defined above tend to increase the energetic (enthalpic) separation between the ground-state and near ground-state conformations on one hand and the open unfolded conformations on the other, or disfavor conformations with intermediate energy (enthalpy), or both. Both of these effects would lead to higher calorimetric cooperativity.4,5 The thermodynamic ramifications of these interactions are further explored in the lower panel of Figure 2, which covers a broad range of values for $`_{\mathrm{coop}}`$ and $`E_{\mathrm{gs}}`$. An additional scenario in which an extra favorable energy for the ground-state is augmented to model (i), viz.,

$$E=E_{\mathrm{contact}}^{}+\gamma _{\mathrm{lh}}N_{\mathrm{lh}}+E_{\mathrm{gs}}$$
$`(4)`$
is also studied \[curve (a)\]. The trend observed in the lower panel of Figure 2 is that stronger many-body cooperative interactions of the type defined above generally lead to higher calorimetric cooperativity. However, there appears to be an upper limit on $`\kappa _2`$ ($`0.96`$) achievable by the helix-packing term alone \[curve (b)\], because at very high $`_{\mathrm{coop}}`$ values it is possible that some intermediate non-ground-state conformations can become relatively stable (c.f. Figure 1).

Figure 3 presents the chevron plots for the three models considered in the upper panel of Figure 2. To model folding and unfolding kinetics at different interaction strengths, an energetic scaling parameter $`ϵ`$ is introduced. At a given $`ϵ`$, the effective energy of a conformation with energy $`E`$ (given by equations 1, 2 or 3) is equal to $`ϵE`$; and variation in $`ϵ/k_\mathrm{B}T`$ (at constant $`T`$) serves as a model denaturant concentration variation, as in ref. 7. Figure 3 shows that at sufficiently strong intrachain interaction (more negative $`ϵ/k_\mathrm{B}T`$), every folding arm of the three chevron plots exhibits a rollover. This suggests that chevron rollover is practically unavoidable in polymer models with physically plausible interactions, because when intrachain interactions become generally very favorable, kinetic trapping is bound to increase in importance.7,8 However, native thermodynamic stability would be extremely high when the model parameter $`ϵ/k_\mathrm{B}T`$ becomes extremely negative. Many such situations are not physically realizable in real proteins,7 whose native stabilities even in zero denaturant are often marginal. In this light, in comparing Figure 3 with experiments, the relevant question is whether there is a quasi-linear regime of the model chevron plots that is consistent with the two-state thermodynamics of the given model and covers a range of thermodynamic stability similar to that of real, simple two-state proteins.

Pursuing this logic, we note that folding rollover occurs quite near to the transition midpoint for model (i), but the $`ϵ/k_\mathrm{B}T`$ range of a quasi-linear regime is more extended for the two more cooperative models (ii) and (iii). The folding arms of models (ii) and (iii) are identical because, by construction, while the $`E_{\mathrm{gs}}`$ term in equation 3 slows unfolding, it does not affect the kinetics of folding. For model (iii), we have used standard histogram techniques and extensive conformational sampling7 at $`ϵ/k_\mathrm{B}T=2.105`$ to determine the dependence of the free energy of unfolding $`\mathrm{\Delta }G_\mathrm{u}`$ on $`ϵ/k_\mathrm{B}T`$ (detailed data not shown). The resulting thermodynamic relation, which is essentially linear and has approximately the same transition midpoint as determined from the kinetic chevron plot, was applied to construct the dotted V-shape in Figure 3. The close agreement between the dotted V-shape and the simulated chevron plot for model (iii) from $`ϵ/k_\mathrm{B}T2.3`$ to $`1.6`$ implies that the folding/unfolding kinetics of model (iii) is consistent with a simple two-state description within this regime. The strongest intrachain interaction in the two-state regime is at $`ϵ/k_\mathrm{B}T2.31`$, which corresponds to a native stability $`\mathrm{\Delta }G_\mathrm{u}10k_\mathrm{B}T`$ for this particular model. It is clear from comparing the chevron plots of models (ii) and (iii) that the linear regime can be readily extended by increasing the magnitude of $`E_{\mathrm{gs}}`$ beyond that in model (iii). But even as it stands, model (iii)’s behavior in Figure 3 should provide a semi-quantitative rationalization for the simple two-state kinetics of many small single-domain proteins. Indeed, more than half of the 24 two-state proteins listed by Plaxco et al. (2000)2 have native stabilities around 25C comparable to or lower than $`10k_\mathrm{B}T`$. For example, $`\mathrm{\Delta }G_\mathrm{u}=3.6k_\mathrm{B}T`$ (2.1 kcal/mol) for CspB at 25C and pH 7.0 (ref. 23), and $`\mathrm{\Delta }G_\mathrm{u}=9.0k_\mathrm{B}T`$ (5.3 kcal/mol) for protein L at 22C and pH 7.0 (ref. 24). In our view, therefore, simple two-state folding/unfolding kinetics emerges as a limiting-case phenomenon when the hypothetically high native stability at which chevron rollover would occur is not attainable by a small single-domain protein. Conversely, rollover becomes observable when a protein fails to achieve a sufficiently high thermodynamic cooperativity commensurate with its native stability.

A comparison between the kinetic properties in Figure 3 and the thermodynamic properties in Figure 2 indicates that the thermodynamic requirement of simple two-state behavior is stringent, allowing only for a small adjustment from empirical baseline subtraction.5 Apparently, a model has to be nearly or as cooperative as model (iii) or more to achieve a reasonable reproduction of simple two-state protein folding/unfolding kinetics. This suggests strongly that, in modeling situations when the heat capacity contributions from bond vibrations are not considered (as in the present cases) and intrachain interaction energies are taken to be temperature-independent, a without-baseline-subtraction $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratio of $`\kappa _2>0.9`$ would likely be required for simple two-state kinetics \[Figure 2, upper panel, model (iii)\].

We have also checked that folding kinetics is essentially single exponential within the quasi-linear regime by confirming that the logarithmic distribution of folding first passage times is approximately linear.7,8,25 In fact, extensive testing for ten values of $`ϵ/k_\mathrm{B}T`$ between $`2.22`$ and $`2.78`$ covering the quasi-linear regime and beyond indicates that they are consistent with single exponential relaxation. Unfolding kinetics is essentially single-exponential as well (detailed data not shown). As in our previous study,7 the onset of non-exponential folding relaxation at interaction strength $`ϵ/k_\mathrm{B}T2.9`$ is concomitant to that of a drastic chevron rollover.

Not surprisingly, Figure 4a shows that the free energy barrier separating the native and denatured states is higher for a more cooperative model, consistent with its slower folding and unfolding rates at the transition midpoint (c.f. Figure 3). Figure 4b shows that the relation between energy and the number of native contacts are approximately linear for the two cooperative models. In this regard, the present exercise suggests that certain many-body interactions embodying a local-nonlocal cooperative interplay ($`_{\mathrm{coop}}=1.0`$) and an added ground-state stability ($`E_{\mathrm{gs}}=2.0`$) in proteins can lead to remarkable improvements in kinetic cooperativity even when the magnitudes of these terms are relatively small.

MODELING A PARTIAL SEPARATION BETWEEN THE
INTERACTIONS FOR THERMODYNAMIC STABILITY AND THE
DRIVING FORCES FOR FOLDING KINETICS

Having established a plausible scenario for simple two-state protein folding/unfolding kinetics, we proceed to broaden our exploration and to better delineate how various energetic components might contribute to this remarkable behavior. As a first step, we consider in this section a somewhat different class of models in which a local-nonlocal cooperative interplay is absent but the unique ground-state conformation is favored by an extra strong energy. The interaction scheme is a simplified version of equation 4 above, with the energy function

$$E=E_{\mathrm{G}\overline{\mathrm{o}}}+E_{\mathrm{gs}},$$
$`(5)`$
where $`E_{\mathrm{G}\overline{\mathrm{o}}}`$ is the usual lattice Gō potential that assigns a favorable energy $`ϵ_0`$ ($`<0`$) to every contact in the native structure of the model and assigns zero energy to all other (nonnative) contacts, and $`E_{\mathrm{gs}}`$ applies only to the ground-state conformation, as in equations 3 and 4. To explore the effect of native topology, we study three cubic-lattice 27mer models with $`E`$ given by equation 5 for three different ground-state structures (Figures 5 and 6).

As discussed above, the $`E_{\mathrm{gs}}`$ term serves to increase native stability and enhance thermodynamic stability, leading to a reduced unfolding rate. But it has no effect on the folding kinetics modeled by Monte Carlo dynamics with the Metropolis acceptance criterion.26,27 Thus, the energetics described by equation 5 entails a partial separation between the interactions that drive the protein to fold kinetically (the pairwise contact energies $`E_{\mathrm{G}\overline{\mathrm{o}}}`$) and the interactions that stabilize the ground-state structure ($`E_{\mathrm{G}\overline{\mathrm{o}}}`$ and the many-body $`E_{\mathrm{gs}}`$). Since $`E_{\mathrm{G}\overline{\mathrm{o}}}`$ contributes partially to native stability, the role separation just described becomes a more predominant feature of the model when $`E_{\mathrm{gs}}`$ is large compare to $`ϵ_0`$. A similar mechanism of partial separation between folding-kinetics and native-stabilizing interactions is also presumed by equation 3 \[model (iii)\] and equation 4 above. Our interest in this scenario was partly motivated by experimental studies showing that mutants of a wildtype protein are much less likely to have a slower unfolding rate than to have a faster folding rate. For example, among the 41 mutants of Fyn SH3 domain studied by Northey et al.,28 only 3 have slightly reduced unfolding rates relative to that of the wildtype, whereas five times as many (15) mutants have folding rates faster than that of the wildtype. This means that interactions that accelerate folding do not necessarily lead to higher native stability (only 4 mutants are more stable than the wildtype), presumably because some mutants that fold fast do not pack well when folded28 (c.f. discussion of conformational strain by Ventura et al.29). These observations suggest that a partial separation of folding-kinetics and native-stabilizing intraprotein interactions as envisioned by the $`E_{\mathrm{gs}}`$ term is physically plausible.

Figure 5 shows that combining a pairwise Gō potential with an $`E_{\mathrm{gs}}`$ term can also lead to simple two-state protein folding/unfolding kinetics, although for these relatively short chains the $`E_{\mathrm{gs}}/ϵ_0`$ ratio needed to achieve simple two-state behavior is large. Figure 5a provides a series of unfolding chevron arms for different $`E_{\mathrm{gs}}`$ values, showing clearly that the quasi-linear regime of the chevron plot can be extended by a more negative $`E_{\mathrm{gs}}`$. For all three models considered in Figure 5 with $`E_{\mathrm{gs}}/ϵ_0=14`$, approximately simple two-state behavior persists to $`ϵ/k_\mathrm{B}T1.15`$ (c.f. simulated chevron plots and dotted V-shapes), corresponding to native stabilities $`\mathrm{\Delta }G_\mathrm{u}10k_\mathrm{B}T`$. As for model (iii) above (Figure 3), for each model in Figure 5, we have verified that folding relaxation is essentially single exponential for $`ϵ/k_\mathrm{B}T>1.6`$ by obtaining linear logarithmic first passage time distributions for six $`ϵ/k_\mathrm{B}T`$ values from $`1.67`$ to $`0.91`$. Unfolding relaxation is also essentially single-exponential (detailed data not shown). Folding relaxation becomes non-exponential (for $`ϵ/k_\mathrm{B}T<1.8`$) only when native stability is much higher than that spanned by the simple two-state regime between $`ϵ/k_\mathrm{B}T1.15`$ and $`ϵ/k_\mathrm{B}T0.723`$.

Figure 6 compares chevron plots of the three cooperative 27mer models. The rank ordering of their folding rates is consistent with a correlation between slower folding rate and higher relative contact order (CO).17 However, for these models, the dependence of folding rate on CO is weak. Near the onset of drastic chevron rollover and non-exponential folding relaxation ($`ϵ/k_\mathrm{B}T1.75`$), the 27mer model with CO $`=0.51`$ folds only approximately 4 times slower than the 27mer model with CO $`=0.28`$. The dispersion in folding rate is even smaller within the simple two-state regime. This is a far cry from the six orders of magnitude of variation in folding rates observed among real, small, single-domain proteins.2 Recently, CO-dependent folding rates have been addressed using explicit-chain models with limited yet encouraging successes. Using a Gō-like potential for 18 small single-domain proteins, Koga and Takada31 obtained a correlation between CO and folding rate, but the variation in rates covered only $`1.5`$ orders of magnitude. More recently, Jewett et al.30 conducted an extensive lattice 27mer simulation study using cooperative Gō-like models with a nonlinear $`E`$$`Q`$ relation. A correlation between CO and folding rate was found but again the dispersion in folding rates spanned only $`1`$ to $`1.5`$ order of magnitude. While the mechanisms and energetics of CO-dependent folding remain to be better elucidated,32 our more recent investigation shows that models with a local-nonlocal cooperative interplay similar to that in models (i) and (ii) (equations 2 and 3) above can lead to a relatively large dispersion in folding rates and a better correlation between CO and folding rate.33,34

A RATIONALIZATION OF NON-ARRHENIUS PROTEIN
FOLDING/UNFOLDING KINETICS

The physics embodied by the extra favorable ground-state energy $`E_{\mathrm{gs}}`$ in the models described by equations 3–5 above implies that there is a fundamental asymmetry between folding and unfolding kinetics.28 This led us to ask whether the same physical picture may shed light on the significant difference in the degree of deviation from Arrhenius kinetics for folding versus unfolding that are often observed in experiments. Early measurements by Segawa and Sugihara35 showed that the folding kinetics of hen egg-white lysozyme was significantly non-Arrhenius (logarithmic folding rate $`\mathrm{ln}k_\mathrm{f}`$ nonlinear in $`1/T`$) whereas the unfolding kinetics was essentially Arrhenius (logarithmic unfolding rate $`\mathrm{ln}k_\mathrm{u}`$ linear in $`1/T`$). Table I summarizes more recent experimental data from the literature for several proteins with simple two-state folding/unfolding kinetics and whose temperature-dependent rates of both folding and unfolding have been measured directly. For the proteins listed, the trend that folding is more non-Arrhenius than unfolding is quantified by reported activation heat capacities for folding $`(\mathrm{\Delta }C_p^{})_\mathrm{f}`$ that are significantly larger in magnitude than the corresponding activation heat capacity for unfolding $`(\mathrm{\Delta }C_p^{})_\mathrm{u}`$. Table I puts the “$`(\mathrm{\Delta }C_p^{})_\mathrm{f}/(\mathrm{\Delta }C_p^{})_\mathrm{u}`$” ratio in quotation marks because the common approach of using temperature-independent activation heat capacities to analyze folding/unfolding kinetics data may be problematic.17 We note that another potential source of the difficulty is that possible temperature dependencies of the heat capacities39 associated with protein folding/unfolding transitions were not considered in such analyses. Nonetheless, as an empirical parameter, “$`(\mathrm{\Delta }C_p^{})_\mathrm{f}/(\mathrm{\Delta }C_p^{})_\mathrm{u}`$” serves well to demonstrate that $`\mathrm{ln}k_\mathrm{f}`$ is often significantly more curvilinear in $`1/T`$ than $`\mathrm{ln}k_\mathrm{u}`$.

This trend can be captured qualitatively in the present modeling context if the solvent-mediated driving forces for folding kinetics and the many-body native-stabilizing interactions are taken to have different temperature dependencies. This is a physically plausible assumption because some intraprotein solvent-mediated forces such as the hydrophobic effect are known to be sensitive to the sizes and shapes of the interacting groups.40-42 Here we use a collection of cooperative 27mer models in Figure 5a with different values of $`E_{\mathrm{gs}}`$ to expound the principles involved. Temperature-dependent interactions are now introduced by letting the pairwise contact energy $`ϵ_0`$ to vary with temperature in a hydrophobic-like manner while leaving $`E_{\mathrm{gs}}`$ temperature-independent.

The following schematic analysis of $`T`$-dependent folding and unfolding rates (Figure 7) is similar to that introduced by Chan26 and Chan and Dill.27 However, the present focus on folding-unfolding asymmetry in the context of three-dimensional protein chain models was not addressed in these earlier studies of short-chain two-dimensional models.26,27 The first step in the present analysis is to obtain from Figure 5a the logarithmic folding and unfolding rates $`\mathrm{ln}k_\mathrm{f}`$ and $`\mathrm{ln}k_\mathrm{u}`$ \[which are taken to be their respective $`\mathrm{ln}(\mathrm{MFPT})`$\] as functions of $`ϵ_0`$ and $`E_{\mathrm{gs}}`$. Since the effective energy is given by $`ϵE`$ in Figure 5 (see above) where $`E`$ is given by equation 5 with $`ϵ_0`$ set to $`1`$, each $`(ϵ,E_{\mathrm{gs}})`$-dependent datapoint in Figure 5a may be regarded as the folding or unfolding rate for the energy function $`E`$ itself with an $`ϵ_0`$ value equals to that of $`ϵ`$ and an $`E_{\mathrm{gs}}`$ value equals $`ϵ`$ times the $`E_{\mathrm{gs}}`$ for the given unfolding chevron arm. We note that within the quasi-linear single-exponential regime,

$$\mathrm{ln}k_\mathrm{f}=\alpha _\mathrm{f}+\beta _\mathrm{f}\frac{ϵ_0}{k_\mathrm{B}T}$$
$`(6)`$
holds approximately for constant $`\alpha _\mathrm{f}`$ and $`\beta _\mathrm{f}`$, because folding kinetics is independent of $`E_{\mathrm{gs}}`$. A least-square fit yields $`\beta _\mathrm{f}=15.4`$. For unfolding within the quasi-linear single-exponential regime, the approximate linear relation

$$\mathrm{ln}k_\mathrm{u}=\alpha _\mathrm{u}+\beta _\mathrm{u}\frac{ϵ_0}{k_\mathrm{B}T}+\beta _\mathrm{u}^{}\frac{E_{\mathrm{gs}}}{k_\mathrm{B}T}$$
$`(7)`$
is expected, where $`\alpha _\mathrm{u}`$, $`\beta _\mathrm{u}`$, and $`\beta _\mathrm{u}^{}`$ are constants. Extensive analyses indicate that $`\beta _\mathrm{u}4.0`$ and $`\beta _\mathrm{u}^{}1.0`$. We use $`\beta _\mathrm{u}=3.9`$, $`\beta _\mathrm{u}^{}=1`$ below. Figure 7a shows that these values fit the simulated unfolding rates of the more cooperative models extremely well.

Next, an hypothetical temperature-dependent $`ϵ_0`$ $`=ϵ_0(T)`$ is introduced in Figure 7b. This functional form for $`ϵ_0/k_\mathrm{B}T`$ (solid curve, left scale) was motivated by the temperature dependence of hydrophobic effects39 and is similar to that explored in refs. 26 and 27. It follows that the temperature-dependent folding rate in the quasi-linear single-exponential regime is now given by equation 6 above with $`ϵ_0`$ $``$ $`ϵ_0(T)`$ provided by Figure 7b, viz.,

$$\mathrm{ln}k_\mathrm{f}(T)=\alpha _\mathrm{f}+\beta _\mathrm{f}\frac{ϵ_0(T)}{k_\mathrm{B}T}.$$
$`(8)`$
Similarly, the temperature-dependent unfolding rate in the quasi-linear single-exponential regime of an effective $`E_{\mathrm{gs}}=14`$ cooperative model is given by equation 7 above with $`ϵ_0`$ $``$ $`ϵ_0(T)`$ provided by Figure 7b while $`E_{\mathrm{gs}}`$ remains temperature independent:

$$\mathrm{ln}k_\mathrm{u}(T)=\alpha _\mathrm{u}+\beta _\mathrm{u}\frac{ϵ_0(T)}{k_\mathrm{B}T}+\frac{E_{\mathrm{gs}}}{k_\mathrm{B}T},$$
$`(9)`$
where the $`E_{\mathrm{gs}}/k_\mathrm{B}T`$ term is set equal to $`14`$ for a reference temperature ($`T^{}`$) in Figure 7c at which $`ϵ_0/k_\mathrm{B}T=1`$. Hence $`E_{\mathrm{gs}}/k_\mathrm{B}T=14(T^{}/T)`$ is linear in $`1/T`$. These temperature-dependent folding and unfolding rates are plotted in the upper part of Figure 7c. It is clear that folding is significantly more non-Arrhenius than unfolding because the only source of non-Arrhenius behavior in the present formulation of the model is $`ϵ_0(T)`$, and $`\mathrm{ln}k_\mathrm{f}`$ depends more strongly on $`ϵ_0(T)`$ ($`\beta _\mathrm{f}=15.4`$) than $`\mathrm{ln}k_\mathrm{u}`$ ($`\beta _\mathrm{u}=3.9`$).

One missing physical ingredient in the consideration thus far is that intrinsic conformational transition rates should accelerate at higher temperature. This is not taken into account if physical time is simply identified with number of attempted Monte Carlo moves, as in the analysis above. The issue has been identified and addressed in some detail in refs. 26 and 27. As in these references, Figure 7b introduces an adjustment factor (dotted line, right scale) to better mimic physical time. Here $`A(T)`$ is the temperature-dependent time needed for a given kinetic process and $`A_0`$ is a reference time. Thus the hypothetical $`\mathrm{ln}[A(T)/A_0]`$ function in Figure 7b stipulates that the intrinsic logarithmic rate \[$`\mathrm{ln}A(T)`$\] is higher at higher temperatures (i.e., the Monte Carlo clock should run faster at higher $`T`$). This adjustment factor is readily incorporated26,27 by setting

| $`\mathrm{ln}(\mathrm{folding}\mathrm{rate})`$ | $`=\mathrm{ln}k_\mathrm{f}(T)\mathrm{ln}[A(T)/A_0],`$ |
| --- | --- |
| $`\mathrm{ln}(\mathrm{unfolding}\mathrm{rate})`$ | $`=\mathrm{ln}k_\mathrm{u}(T)\mathrm{ln}[A(T)/A_0],`$ |

$`(10)`$
where $`\mathrm{ln}k_\mathrm{f}`$ and $`\mathrm{ln}k_\mathrm{u}`$ on the right hand side are the expressions given respectively by equations 8 and 9. Macroscopically, this amounts to introducing an additional enthalpic contribution24,26 to the free energy barrier of protein folding. The lower part of Figure 7c shows that incorporating a $`\mathrm{ln}[A(T)/A_0]`$ term can lead to a more realistic description of temperature-dependent protein folding and unfolding rates (c.f. experimental data in Figure 7d; see also Figure 3 of ref. 23, Figure 1C of ref. 24, Figure 3 of ref. 36, and Figure 5 of ref. 38).

DISCUSSION: A NEAR-LEVINTHAL SCENARIO FOR SIMPLE
TWO-STATE PROTEIN FOLDING/UNFOLDING KINETICS

The results of the present investigation suggest strongly that the physical interactions underlying the simple two-state folding/unfolding kinetics of small single-domain proteins should involve many-body effects beyond that stipulated by common additive Gō models, even though the physico-chemical origins of these effects remain to be elucidated. Therefore, with regard to protein thermodynamic and kinetic cooperativity, common Gō models with pairwise additive contact energies are not ideal. Apparently, the type of many-body interactions that are conducive to simple two-state kinetics also lead to higher thermodynamic cooperativity, and entail a partial separation between folding-kinetics and native-stabilizing interactions. Historically, an impetus to formulate the Levinthal paradox might have been the discovery in the late 1960s that some proteins were calorimetrically two-state43 (see ref. 4 and references therein). Naturally, an extreme interpretation of the $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}=1`$ property would imply that only two enthalpy levels exist (native and denatured), and thus the landscape should resemble a golf course (Figure 8a). However, a golf-course landscape dictates that folding would be exceedingly slow, but the folding of real proteins is relatively fast. To address the “why is folding fast” question, recent theoretical discussions emphasize the funnel-like nature of the protein folding energy landscape as a solution to the Levinthal paradox44,45 (Figure 8b); and common Gō potentials are often used to model a relatively smooth funnel-like energy landscape.12 We found that calorimetric two-state cooperativity can be consistent with the funnel-like landscapes of three-dimensional Gō models, provided some lattitude is allowed for empirical baseline subtractions. This is because in these native-centric models, the conformational populations with intermediate energies (enthalpies) — though not zero — are relatively low.5-8

However, as discussed above, common Gō models are insufficient for simple two-state protein folding kinetics.7,8 Because kinetic traps are still significant in these constructs under native conditions, the chevron plots they predicted have severe rollovers.46 Another shortcoming of common Gō models is that their predicted folding rates are often too fast compared to that of real proteins.8,14 In one example, it is at least four orders of magnitude faster.8 So, in a sense, in the context of recent Gō modeling efforts, the critical question has shifted from “why is folding fast” to “why is folding slow.” The present study concludes that a thermodynamic cooperativity higher than that afforded by common additive Gō models is necessary for simple two-state kinetics (Figures 2 and 3). This scenario also offers a clue to the “why is folding slow” question. For real, small, single-domain proteins, it appears that the key to avoiding kinetic traps and chevron rollover is to have only weakly favorable intrachain interactions during the folding process (gentle upper slopes of the funnel in Figure 8c) until a significant fraction of the chain is native-like and ready to come together to form a large number of native contacts at once, at which point strong cooperative many-body effects kick in to stabilize the structure (steep lower slopes of the funnel in Figure 8c). This idea is implemented in the cooperative models studied here. Indeed, in Figures 3 and 5, the quasi-linear folding regimes of the cooperative models are in the weakly-interacting unfolding regimes (small $`ϵ/k_\mathrm{B}T`$) of the corresponding additive Gō models. Thus, for folding of the cooperative models, the energetic bias is not very strong during most of the conformational search. This feature serves to diminish the effects of kinetic traps (shallow minima on the gentle upper slopes of Figure 8c), because the depths of kinetic traps in heteropolymers are often correlated with the overall energetic bias towards the native structure.27 As a result, folding is faster in the cooperative models relative to other heteropolymer models with deep kinetic traps. But at the same time, the very feature of a weakened energetic bias towards the native structure during most of the conformational search also leads to slower folding in comparison with common Gō models, because the latter have stronger native biases during the corresponding kinetic process. Nonetheless, the reduction in folding rate relative to common Gō models does not make the cooperative models less proteinlike, because even a small bias is sufficient to circumvent the Levinthal paradox, as all of our models fold. In this way, the cooperative model scenario provides a physical plausible answer to the “why is folding slow” question posed above. In fact, for real proteins, folding rates may be further reduced by the inevitable presence of favorable nonnative interactions, which are not taken into account by the native-centric cooperative models here. Anti-cooperativity of certain hydrophobic interactions47 may also play a role in discouraging premature chain collapse. In this scenario, the energy landscape of a simple two-state protein is still funnel-like but with a narrow bottleneck (Figure 8c). The resulting highly bimodal distribution of energy and high thermodynamic cooperativity thus approach (though never equal) that of an hypothetical Levinthal golf course. We emphasize that the relatively smooth funnel drawings in Figure 8 should be viewed only as pictorial devices for underscoring the smoothness of the energy landscapes of native-centric models relative to that of models with deeper kinetic traps. Even for Gō and Gō-like models, energy landscapes cannot be completely smooth because of repulsive interactions (including excluded volume effects) and other microscopic energy barriers due to bond rotations and solvation, for example (refs. 8, 26, 27; c.f. equation 10).

Part of the present view is similar to that of Jewett et al.,30 who recently introduced a native-centric model in which the energy $`E`$ of a conformation — unlike that in common Gō models — does not decrease linearly with the number of native contacts. (Conformations with more negative $`E`$’s are more favorable.) In their model, the rate of decrease in $`E`$ with increasing fractional number of native contact $`Q`$ becomes progressively higher as the native structure is approached (when $`Q1`$). This means that the energetic bias towards the native structure is not strong during the initial stages of folding when there are few native contacts (small $`Q`$), but becomes stronger when the native structure is approached during the final stages of folding ($`Q1`$). Insofar as this general trend is concerned, the physical picture of cooperative folding discussed above (especially the model defined by equation 5) is very much similar to that of Jewett et al. Nevertheless, although both the model of Jewett et al. and the present cooperative models have high degrees of thermodynamic cooperativity, their underlying mechanisms are not identical. More recent investigations indicate that detailed kinetic features, such as the correlation between CO and folding rate, do depend significantly on how thermodynamic cooperativity is achieved microscopically.33,34 Since the lattice model of Jewett et al. was inspired by the more general topomer search model of folding,48 it would be extremely interesting to compare in future investigations the relationships between the topomer search model and the several different scenarios of thermodynamic cooperativity explored here.

As we have emphasized,4,5 the experimental calorimetric two-state criterion, which has proven useful for evaluating protein chain models,4-8,33,34,49-52 does not imply that there are only two infinitely sharp energy (or enthalpy) levels. In other words, the calorimetric two-state criterion does not preclude the existence of “partially unfolded” conformations with energies intermediate between the energy distribution peaks under strongly folding and strongly denaturing conditions (c.f. Figure 16 of ref. 4, Figure 10 of ref. 5, Figure 1 of ref. 7, and Figures 7–9 of ref. 8). These theoretical findings are consistent with native-state hydrogen exchange experiments.53 Conformational populations with intermediate energies in several calorimetrically cooperative models tested thus far (see figure references above) share some similarities with that in calorimetrically non-cooperative constructs such as certain HP models4 and a 15mer 20-letter sidechain model.5,54,55 However, the critical difference between calorimetric cooperative and non-cooperative models is that, at the transition midpoint, conformational populations with intermediate energies are not significant for calorimetrically cooperative models but are significant for calorimetrically non-cooperative models. This difference is well characterized by the $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratio.4,5,49-52 Our approach of evaluating chain models by experimental cooperativity criteria is designed to address what kind of elementary intrachain interactions may be needed to produce the generic cooperative features of proteins, while taking into account as much as possible that proteins are polymers and therefore chain connectivity and excluded volume are severe constraints. In this respect, investigations using self-contained polymer models such as that conducted here are fundamentally more informative than those that use postulated free energy profiles in the absence of explicit chain representations (see, e.g., ref. 56).

In summary, our results suggest that intramolecular recognition in real two-state proteins is highly specific. As well, the role of many-body interactions in providing a larger average energetic difference between “native” and “denatured” conformations than that afforded by common pairwise additive interaction schemes have potentially important implications for the discrimination of decoys in protein structure prediction.57 In principle, the many-body interactions proposed in the present work can be characterized quantitatively through careful experiments and extensive atomic simulations. How sidechain packing, sidechain/mainchain correlation,58 and interactions such as hydrogen bonding contribute to this mechanism remains to be investigated. To help address these questions, the ramifications of the different scenarios explored in this work need to be first delineated in greater detail.33,34 Lattice model studies of Klimov and Thirumalai showed that sidechain degrees of freedom increase the sharpness of the thermodynamic folding/unfolding transition relative to that of the corresponding (mainchain) model with no sidechains.54 However, their short-chain 20-letter sidechain models configured on three-dimensional simple cubic lattices do not appear to be calorimetrically cooperative. The $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratio of one of the sidechain sequences studied in refs. 54 and 55 is equal to $`\kappa _2=0.38`$ without baseline subtraction, and increases only to $`\kappa _2^{(\mathrm{s})}=0.54`$ after reasonable baseline subtractions.5 These values are far from the $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}1`$ required for calorimetrically cooperative behavior. Kinetically, even a Gō-like version of their sidechain model exhibits a severe chevron rollover (Figure 3 of ref. 59). These results imply that while sidechain contributions are expected to be important for protein cooperativities,54 their role has yet to be better elucidated.

Finally, while the present study focuses on the behavior of small single-domain proteins, we hasten to emphasize that not all proteins have simple two-state folding/unfolding kinetics. Hence, the high cooperativity requirement deduced in the above analysis may not apply to other proteins. In fact, one distinct advantage of the energy landscape perspective and self-contained polymer modeling is their ability to cover a wide spectrum of possible protein behavior under a unified physical framework. For instance, although common additive Gō models are insufficient for simple two-state kinetics, they are useful for understanding real proteins with chevron rollovers.7,8,46,60 And even calorimetrically non-cooperative models (see discussion in ref. 5) may prove to be helpful in rationalizing downhill protein folding61 as well.

Acknowledgments

We thank Yawen Bai, Robert L. Baldwin, Alan Davidson, Julie Forman-Kay, Michael Levitt, Kevin Plaxco, Steve Plotkin, Boris Steipe, and Dev Thirumalai for helpful discussions, and Kevin Plaxco for kindly providing ref. 30 before publication. This work was partially supported by the Canadian Institutes of Health Research (CIHR grant no. MOP-15323), a Premier’s Research Excellence Award from the Province of Ontario, and the Ontario Centre for Genomic Computing at the Hospital for Sick Children in Toronto. H. S. C. is a Canada Research Chair in Biochemistry.

References

Table I. Deviations from Arrhenius behavior in protein folding and unfolding.

| Protein name | $`(\mathrm{\Delta }C_p^{})_\mathrm{f}/(\mathrm{\Delta }C_p^{})_\mathrm{u}`$ |
| --- | --- |
| T4 lysozyme mutanta | $`2.99`$ |
| CI2b | $`2.95`$ |
| CspBc | $`9.0`$ |
| Protein Ld | $`1.68`$ |
| NTL9e | $`1.33`$ |

a Chen et al. (ref. 36)
b Jackson & Fersht (ref. 37)
c Schindler & Schmid (ref. 23)
d Scalley & Baker (ref. 24)
e Kuhlman et al. (ref. 38)

Figure Captions

Figure 1. Modeling a cooperative interplay between local conformational preference and protein core packing in a 55mer native-centric four-helix-bundle model. Certain native helices are shown as thick solid lines for illustrative purposes. Their sequence positions are indicated by the thin dotted lines depicting the rest of the full native structure. Given a native helix is completely formed \[the front right helix appearing in both (a) and (b) in this example\], a favorable cooperative energy $`_{\mathrm{coop}}`$ ($`<0`$) is assigned if either (a) 4 or more of the 6 native contacts (thick dotted lines) between the given fully formed native helix and each of the two chain segments for the two flanking native helices are present, or (b) at least 2 of the 3 “diagonal” native contacts are present between the given fully formed native helix and the chain segment for the diagonally neighboring native helix, or both \[(a) and (b)\]. Thus, condition (a) requires at least 8 inter-helix nearest-neighbor native contacts, whereas condition (b) requires at least 2 inter-helix next-nearest-neighbor native contacts on the simple cubic lattice. It follows that the maximum total contribution from these cooperative energies is $`4_{\mathrm{coop}}`$ when all four native helices are completely formed and correctly packed against one another.

Figure 2. The overall thermodynamic cooperativity of a model protein is boosted by the many-body interactions described in Figure 1 and in the text. Upper panel: Heat capacity as a function of temperature for (i) the pairwise additive native-centric 55mer model (equation 1), (ii) a model having the interactions in (i) plus the cooperative interaction scheme in Figure 1 with $`_{\mathrm{coop}}=1.0`$ (equation 2), and (iii) a model with the interactions in (ii) plus an extra favorable energy of $`E_{\mathrm{gs}}=2.0`$ for the ground-state conformation (equation 3). The heat capacity scans were obtained using Monte Carlo histogram techniques based on conformational sampling around each model’s transition midpoint.5-7 The inset (from ref. 7) shows the 55mer model ground-state structure with (nominally) hydrophobic and polar residues depicted respectively as filled and open circles. Lower panel: Van’t Hoff to calorimetric enthalpy $`\mathrm{\Delta }H_{\mathrm{vH}}/\mathrm{\Delta }H_{\mathrm{cal}}`$ ratios are given by $`\kappa _2`$ defined in ref. 5 (without empirical baseline subtractions) for three classes of 55mer models whose interaction schemes are parametrized by an $`_{\mathrm{coop}}`$ variable. (a) As in (i) above plus an extra favorable energy of $`_{\mathrm{coop}}`$ for the ground-state conformation (equation 4). (b) As in (ii) above but, instead of fixing $`_{\mathrm{coop}}=1.0`$, a variable $`_{\mathrm{coop}}`$ is used for the helix packing contribution defined by Figure 1 and equation 2. (c) As in (b) plus an extra favorable energy of $`E_{\mathrm{gs}}=_{\mathrm{coop}}`$ for the ground-state conformation (equation 3).

Figure 3. Chevron plots for the 55mer models (i), (ii), and (iii) in Figure 2 are provided by the negative natural logarithm of mean first passage time (MFPT) as functions of interaction strength $`ϵ/k_\mathrm{B}T`$. The present Monte Carlo (MC) dynamics simulations use the same general procedure as that in ref. 7, now with a move set consisting of end flips (4.7%), corner flips (58.3%), crankshafts (27%) and rigid rotations (10%). Folding (filled symbols) starts from a randomly generated conformation. First passage time (FPT) for folding is the number of attempted MC moves needed to reach the ground-state conformation. Unfolding (open symbols) starts from the ground-state conformation. Here unfolding FPT is the number of attempted MC moves needed to reach a conformation with fewer than 7 native contacts. Each plotted MFPT is averaged from 1,000 trajectories. Solid and dashed curves through the data points are mere guides for the eye. The dotted V-shape for model (iii) is an hypothetical simple two-state chevron plot consistent with the $`ϵ/k_\mathrm{B}T`$ dependence of thermodynamic stability as given by the free energy difference $`\mathrm{\Delta }G_\mathrm{u}`$ between the ground state and the unfolded conformational ensemble with $`<7`$ native contacts.

Figure 4. Distribution of native contacts in the 55mer models. The total (maximum) number of spatial nearest-neighbor native contacts in the ground-state conformation is 60 (diagonal contacts in Figure 1b are not included in this accounting). $`Q`$ is the fractional number of native contacts for a conformation, defined as the number of native contacts it contains divided by the maximum. (a) The free energy profiles of the models (i), (ii) and (iii) defined in Figure 2 are given by the negative logarithmic distributions of $`Q`$. The profiles shown are for $`ϵ/k_\mathrm{B}T=`$ $`2.34`$, $`2.11`$, and $`2.0`$ respectively for (i), (ii) and (iii), near each model’s transition midpoint, and were obtained by standard MC histogram techniques.5-7 (b) Correlation between energy $`E`$ and $`Q`$ for models (ii) and (iii). Dots indicate the existence of conformations with the given $`E`$ and $`Q`$ values. The open diamond and square mark the ground-state energies of $`40.1`$ and $`42.1`$ for models (ii) and (iii) respectively.

Figure 5. Chevron plots for three different 27mer Gō models and for their corresponding models with an extra favorable energy $`E_{\mathrm{gs}}`$ assigned to the native (ground-state) structure (as shown), with $`ϵ_0=1`$ (equation 5). Folding MFPT is independent of $`E_{\mathrm{gs}}`$, and $`ϵ`$ has the same meaning as in Figure 3. Here each MFPT is averaged from 500 trajectories. Folding (filled symbols) and unfolding (open symbols) simulations were performed as for Figure 3 except only local chain moves were used for the present 27mer MC dynamics simulation (no rigid rotations), and unfolding FPT is now defined by the time needed to reach a conformation with fewer than 4 native contacts from a given starting ground-state conformation with 28 native contacts. Solid curves are mere guides for the eye. (a) Unfolding chevron arms for eight models with different degrees of cooperativity are shown; from top to bottom, $`E_{\mathrm{gs}}=`$ $`0`$ (common additive Gō model), $`E_{\mathrm{gs}}=`$ $`2`$, $`4`$, $`6`$, $`8`$, $`10`$, $`12`$, and $`14`$. (b) and (c) Unfolding chevron arms for the common additive Gō model ($`E_{\mathrm{gs}}=0`$, upper curves) are compared with the unfolding arms for $`E_{\mathrm{gs}}=14`$ (lower quasi-linear curves). The dotted V-shapes are hypothetical simple two-state chevron plots consistent with the $`ϵ/k_\mathrm{B}T`$ dependence of $`\mathrm{\Delta }G_\mathrm{u}`$ between the ground state and the unfolded conformational ensemble with $`<4`$ native contacts of a given model; $`\mathrm{\Delta }G_\mathrm{u}`$ values are determined by standard histogram techniques5-7 based on conformational sampling at $`ϵ/k_\mathrm{B}T=0.935`$, $`0.917`$ and $`0.917`$ for (a), (b) and (c) respectively (detailed data not shown).

Figure 6. Comparing chevron plots for the three different models in Figure 5 with $`E_{\mathrm{gs}}=14`$ shows that a model with lower native contact order (CO, as defined in ref. 17) tends to fold slightly faster. Here CO $`=`$ $`0.28`$, $`0.40`$, and $`0.51`$ for models (a), (b) and (c) respectively.

Figure 7. Rationalizing non-Arrhenius protein folding and unfolding kinetics. (a) A linear fit of the logarithmic unfolding rates (vertical axis) of the 27mer models in Figure 5a with $`E_{\mathrm{gs}}6`$ to the expression shown (horizontal axis). For the data points plotted, the correlation coefficient $`r=0.998`$. (b) Solid curve: an hypothetical hydrophobic-like temperature dependence of the model interaction strength $`ϵ_0/k_\mathrm{B}T`$ that drives folding kinetics (left vertical scale). Dashed line: an hypothetical temperature dependence of the intrinsic conformational transition rate $`1/A(T)`$ relative to a reference rate $`1/A_0`$ (right vertical scale). More analytical details of this physical picture are provided in refs. 26 and 27. (Note that a typographical error should be corrected in the caption for Figure 4 of ref. 26: “$`\mathrm{\Delta }S_0=5.2`$” should read “$`\mathrm{\Delta }S_0=5.2`$.”) (c) Upper curves: Temperature-dependent folding and unfolding rates obtained by combining the solid curve in (b) and data from Figure 5a (equations 8 and 9, temperature-dependent intrinsic conformational transition rate not taken into account). Lower curves: Temperature-dependent folding and unfolding rates obtained by combining the solid curve and dashed line in (b) and data from Figure 5a (equation 10, temperature-dependent intrinsic conformational transition rate taken into account). (d) Included for comparison are the temperature-dependent CI2 folding and unfolding rates at 25C and pH 6.3 from the experiments of Jackson and Fersht (adapted from Figure 4 of ref. 37). The upper and lower folding curves were for 0 M and 1.5 M GdnHCl respectively, the unfolding curve was measured at 0 M GdnHCl.

Figure 8. Schematics of hypothetical and proposed energy landscapes for protein folding. (a) A golf-course or “Levinthal” landscape. (b) A funnel landscape. (c) A “near-Levinthal” scenario.


format: text

text:
# Comment to “Experimental Evidence of a Dynamic Jahn-Teller Effect in C⁺₆₀”

A recent Letter S. E. Canton et al., Phys. Rev. Lett. 89, 045502 () (2002) reports photoemission (PE) data for the free C60 molecule, showing an interesting three-peak structure, presented as evidence of dynamic Jahn-Teller (DJT) effect in the C$`{}_{}{}^{+}{}_{60}{}^{}`$ ion. Those data constitute, along with earlier spectra by the Uppsala group P. Brühwiler et al., Chem. Phys. Lett. 279, 85 () (1997), the best available piece of experimental evidence about the spectrum of a hole in fullerene. DJT must indeed affect the fivefold-degenerate $`h_u`$ hole molecular orbital N. Manini et al., Philos. Mag. B 81, 793 () (2001), but we contend that the energy separation of these three peaks is far too large for the proposed tunneling interpretation to be correct.

In detail, the observed structure is claimed to indicate a $`D_{3d}`$ distortion in the lowest-energy JT well, accompanied by tunneling between the local minima. The main argument offered is that the ten $`D_{3d}`$ valleys yield three tunnel-split states $`H_u`$ (5) + $`G_u`$ (4) + $`A_u`$ (1), to be identified with the three observed structures. Alternatively, six $`D_{5d}`$ valleys would split into two $`H_u`$ (5) + $`A_u`$ (1), while valleys of lower symmetry ($`D_{2h}`$, $`C_{2h}`$) would split more than threefold. Given this interpretation however, the tunnel splitting magnitudes between $`D_{3d}`$ valleys would amount to 230 meV ($`H_uG_u`$), and 390 meV ($`H_uA_u`$). Both splittings are alarmingly larger than any of the vibrational frequencies of fullerene ($`32÷195`$ meV).

Our contention is precisely that these supposed splittings are far too large for the tunneling interpretation to be valid. Tunnel splittings only make sense in the large-barrier limit, when they are smaller than the smallest vibrational quantum $`\mathrm{}\omega `$, here of 32 meV (Fig. 1A), and the model of cited Ref.  applies. In the opposite limit (Fig. 1B), barriers between valleys are lower than the kinetic energy, and tunneling-split levels are replaced by extended anharmonic excitations, delocalized over all valleys. To estimate the two lowest excitations to be expected for C$`{}_{}{}^{+}{}_{60}{}^{}`$, we carried out a realistic calculation of the lowest $`H_u`$, $`G_u`$, and $`A_u`$ vibronic states of an $`h_u`$ electronic level JT coupled to eight fivefold $`H_g`$ and to six fourfold $`G_g`$ vibrational modes N. Manini et al., Philos. Mag. B 81, 793 () (2001). Using the ab-initio JT coupling parameters and frequencies of Ref. N. Manini et al., Philos. Mag. B 81, 793 () (2001) in a symmetry-restricted Lanczos diagonalization, we obtain 18 meV and 30 meV ($`\mathrm{}\omega `$) for the excitation energies from the $`H_u`$ ground state to the lowest $`A_u`$ and $`G_u`$ vibronic states respectively. These values are one order of magnitude smaller than the PE structures S. E. Canton et al., Phys. Rev. Lett. 89, 045502 () (2002).

Is this discrepancy due to uncertainty in the precise values of the coupling parameters? We think not. If the actual couplings were smaller or, as is more likely, larger than the those assumed, or if they made $`D_{3d}`$ wells lower than $`D_{5d}`$ wells, then the lowest $`A_u`$, $`G_u`$, and $`H_u`$ states would always lie within a range of $`\mathrm{}\omega 30`$ meV. The spectral structures above 200 meV might reflect high frequency vibrons, or else they might be due to an electronic splitting, in analogy to the interpretation of the PES of Fe(CO)5 in cited Ref. . The observed splittings would be in the same range ($`180`$ meV) as the computed electronic JT splitting of Fig. 3 of Ref. N. Manini et al., Philos. Mag. B 81, 793 () (2001). A full calculation of the spectrum that will include both low-energy (tunneling) and high-energy “electronic” splittings could in principle be done based on Fermi’s golden rule L. Martinelli et al., Phys. Rev. B 43, 8395 ; G. Grosso and G. Pastori Parravicini, Solid State Physics (Academic Press, San Diego, 2000), Ch. XII.7() (1991), but is at present still unavailable.

In summary, it would seem that the three-peak structure, though certainly related to DJT, does not particularly provide evidence for tunneling among $`D_{3d}`$ valleys as claimed. A full explanation of this spectrum must await further work.


text:
# Quantised Three-Pillar Problem

## Abstract

This paper examines the quantum mechanical system that arises when one quantises a classical mechanical configuration described by an underdetermined system of equations. Specifically, we consider the well-known problem in classical mechanics in which a beam is supported by three identical rigid pillars. For this problem it is not possible to calculate uniquely the forces supplied by each pillar. However, if the pillars are replaced by springs, then the forces are uniquely determined. The three-pillar problem and its associated indeterminacy is recovered in the limit as the spring constant tends to infinity. In this paper the spring version of the problem is quantised as a constrained dynamical system. It is then shown that as the spring constant becomes large, the quantum analog of the ambiguity reemerges as a kind of quantum anomaly.

: J. Phys. A: Math. Gen.
Consider a rigid beam of weight $`W`$ supported by three identical incompressible pillars (see Fig. 1). Let the upward forces provided by the pillars be $`F_1`$, $`F_2`$, and $`F_3`$, respectively. If the beam is in static equilibrium, then the sum of the upward forces must equal $`W`$:

$`F_1+F_2+F_3=W.`$ (1)
Also, the torque on the beam about any point must vanish. If we calculate the torque about the centre, we obtain the condition

$`F_1=F_3.`$ (2)
(Calculating the torque about any other point does not give additional information). The simultaneous solution to (1) and (2) is

$`F_1`$ $`=`$ $`{\displaystyle \frac{1}{2}}(Wf),`$
$`F_2`$ $`=`$ $`f,`$
$`F_3`$ $`=`$ $`{\displaystyle \frac{1}{2}}(Wf),`$ (3)
where $`f`$ is an arbitrary force that cannot be determined by the conditions of the problem. This is an elementary example of a classical mechanical system whose physical characteristics are underdetermined.

It is possible to reformulate the three-pillar problem in such a way as to remove this ambiguity. We replace the three pillars by three identical springs having spring constant $`k`$. Now, the rigid beam rests on these springs, and the springs are displaced from their equilibrium lengths by the amounts $`x`$, $`y`$, and $`z`$ (see Fig. 2). For this problem we can determine $`x`$, $`y`$, and $`z`$, and thus we can determine the forces uniquely. The force balance equation reads

$`kx+ky+kz=W,`$ (4)
and the torque balance condition implies that

$`x=z.`$ (5)
The condition that the beam is straight and rigid implies that

$`x+z=2y,`$ (6)
which is a new equation having no analogue in the three-pillar problem. The simultaneous solution of equations (4) – (6) gives $`x=y=z=\frac{W}{3k}`$, and thus the upward forces imposed on the beam by the three springs are

$`F_1={\displaystyle \frac{1}{3}}W,F_2={\displaystyle \frac{1}{3}}W,F_3={\displaystyle \frac{1}{3}}W.`$ (7)
The ambiguity in the three-pillar problem is eliminated because the flexibility in the springs and the rigidity of the beam give rise to the additional condition (6) that allows us to solve for the forces. Thus, if the pillars are replaced by compressible objects, which need not even be identical, then the indeterminacy is lifted. However, if we take the limit $`k\mathrm{}`$ in the spring problem, then the springs become rigid objects, and we recover the three-pillar problem. To see this, we note from (4) that the limit $`k\mathrm{}`$ gives

$`x+y+z=0`$ (8)
because $`W`$ is a constant. The torque condition (5) is still valid, and along with (6) we find that $`x=y=z=0`$. We obtain this result because in this limit the springs become incompressible, and the deviation from the equilibrium length of the springs must vanish. However, the forces are no longer determined because they are expressed in the ambiguous form $`F_i=\mathrm{}\times 0`$ ($`i=1,2,3`$). Thus, any values for the forces are allowed subject to constraints (1) and (2).

Let us now consider the quantum mechanical version of the three-spring problem. To describe the quantised three-spring problem we start with the Hamiltonian

$`H={\displaystyle \frac{1}{2m}}\left(p_x^2+p_y^2+p_z^2\right)+{\displaystyle \frac{1}{2}}k\left(x^2+y^2+z^2\right),`$ (9)
which represents three uncoupled harmonic oscillators. Note that we have shifted the zeros of the variables $`x`$, $`y`$, and $`z`$ by the amount $`\frac{W}{3k}`$ so that the beam oscillates about its classical resting position. We then impose the constraint (6) to eliminate the dynamical variable $`y`$. The Hamiltonian thus obtained is

$`H={\displaystyle \frac{1}{2m}}\left[p_x^2+{\displaystyle \frac{1}{2}}\left(p_x+p_z\right)^2+p_z^2\right]+{\displaystyle \frac{1}{2}}k\left[x^2+{\displaystyle \frac{1}{4}}(x+z)^2+z^2\right].`$ (10)
It is convenient to make the change of variables $`r=x+z`$ and $`s=xz`$. The variable $`r`$ represents vertical oscillatory motion of the centre of mass of the beam and the variable $`s`$ is associated with the rotational motion of the beam. In terms of these variables, the resulting Hamiltonian is diagonal and we obtain the time-independent Schrödinger equation

$`\left[{\displaystyle \frac{1}{2m}}\left(3{\displaystyle \frac{^2}{r^2}}+2{\displaystyle \frac{^2}{s^2}}\right)+{\displaystyle \frac{1}{8}}k\left(3r^2+2s^2\right)\right]\psi (r,s)=E\psi (r,s).`$ (11)
Let us calculate the expectation value of the force, given by the negative derivative of the potential, applied by the central oscillator. We perform this calculation in the ground state of the quantum system. The (unnormalised) ground-state wave function is

$`\psi _0(r,s)=\mathrm{exp}\left({\displaystyle \frac{1}{4}}\sqrt{km}r^2\right)\mathrm{exp}\left({\displaystyle \frac{1}{4}}\sqrt{km}s^2\right).`$ (12)
The expectation value of force in the central oscillator, taking into account the shift of the zero, is given by $`W/3`$ plus the integral

$`0|\mathrm{Force}|0`$ $`=`$ $`{\displaystyle \frac{k}{2}}{\displaystyle \frac{drrds\mathrm{exp}\left[\frac{1}{2}\sqrt{mk}\left(r^2+s^2\right)\right]}{drds\mathrm{exp}\left[\frac{1}{2}\sqrt{mk}\left(r^2+s^2\right)\right]}}`$ (13)
$`=`$ $`{\displaystyle \frac{k}{2}}{\displaystyle \frac{drr\mathrm{exp}\left(\frac{1}{2}\sqrt{mk}r^2\right)}{dr\mathrm{exp}\left(\frac{1}{2}\sqrt{mk}r^2\right)}},`$
where we have divided out the integral over the variable $`s`$. So long as the spring constant $`k`$ is finite, a reflection symmetry argument implies that the integral in the numerator vanishes. Indeed, in any state, the expectation value of $`r`$ vanishes by reflection symmetry. However, as $`k\mathrm{}`$, the expression in (13) is the product of 0 and $`\mathrm{}`$, and we have obtained an ambiguous result.

This ambiguity may be regarded as a kind of quantum anomaly. Typically, in quantum field theory one encounters anomalies for which there is a divergent integral multiplying a quantity that vanishes because of a geometrical symmetry argument. For example, in the Schwinger model of two-dimensional electrodynamics the trace of the photon propagator formally vanishes by the symmetry properties of two-dimensional gamma matrices . However, the photon propagator contains an integral that is logarithmically divergent. There are many tricks that can be used to regulate and then calculate such ambiguous quantities. In the case of the Schwinger model one can regulate the divergent integral by performing the calculation in $`2+ϵ`$ dimensions; in $`2+ϵ`$ dimensions the trace of the gamma matrices does not vanish and the integral is finite. One then takes the limit $`ϵ0`$ to calculate the value of the anomaly. Another example of an anomaly is the so-called trace anomaly, where the stress-energy tensor formally has a vanishing trace, but is represented by a divergent integral .

We view the expression in (13) as a kind of quantum anomaly because there is one quantity (here, the integral) that vanishes by a geometrical symmetry argument (reflection of $`r`$). This quantity is multiplied by a second factor $`k`$ that diverges (here, $`k\mathrm{}`$). Our objective is to demonstrate that, depending on the regularisation scheme chosen, we can get any result for the expectation value (13) of the force.

As an example, we can regulate the anomalous ambiguous product

$`A=\underset{\mathrm{\Lambda }\mathrm{}}{lim}\mathrm{\Lambda }{\displaystyle _{\mathrm{}}^{\mathrm{}}}drre^{r^2}`$ (14)
as follows:

$`A\underset{\mathrm{\Lambda }\mathrm{}}{lim}\mathrm{\Lambda }{\displaystyle _{\mathrm{}}^{\sqrt{\mathrm{log}[\mathrm{\Lambda }/(2a)]}}}drre^{r^2},`$ (15)
where $`a>0`$ is arbitrary. Evaluating (15) exactly, we obtain

$`A=a.`$ (16)
Apparently, by choosing an appropriate regulation scheme we can obtain any value $`a`$ for $`A`$. Since the expression for the force in (13) has the form of the anomalous product in (14), we see that as the spring constant $`k`$ tends to infinity, the expectation value of the force is arbitrary.

We wish to express our gratitude to H. F. Jones for his stimulating suggestions. DCB gratefully acknowledges financial support from The Royal Society. This work was supported in part by the U.S. Department of Energy.