New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
A statistical mechanical model for βhairpin kinetics
Abstract
Understanding the mechanism of protein secondary structure formation is an essential part of the proteinfolding puzzle. Here, we describe a simple statistical mechanical model for the formation of a βhairpin, the minimal structural element of the antiparallel βpleated sheet. The model accurately describes the thermodynamic and kinetic behavior of a 16residue, βhairpinforming peptide, successfully explaining its twostate behavior and apparent negative activation energy for folding. The model classifies structures according to their backbone conformation, defined by 15 pairs of dihedral angles, and is further simplified by considering only the 120 structures with contiguous stretches of native pairs of backbone dihedral angles. This single sequence approximation is tested by comparison with a more complete model that includes the 2^{15} possible conformations and 15 × 2^{15} possible kinetic transitions. Finally, we use the model to predict the equilibrium unfolding curves and kinetics for several variants of the βhairpin peptide.
As is evident from the presentations at this Colloquium, the continuous discovery of thousands of new gene sequences is producing a revolution in all aspects of protein physics, chemistry, and biology. Foremost among these is the proteinfolding problem. C. B. Anfinsen, in his Nobel Prize winning experiments at the National Institutes of Health (1), showed that a denatured protein can refold spontaneously to form a biologically functional (native) structure. From this result, Anfinsen concluded that the information for determining the threedimensional structure is somehow encoded in the amino acid sequence. This work has led to the realization that it should in principle be possible to calculate the threedimensional structure of a protein from its amino acid sequence. Calculating the structure from the sequence has become known as the first part of the proteinfolding problem and currently engages a large number of theoretical and computational scientists. The second part of the proteinfolding problem is to understand how a protein folds. That is, what are the kinetics and mechanism (or mechanisms) of protein folding? This question is in many ways more challenging because for in vitro folding the ultimate answer is a description of the distribution of threedimensional structures as a function of time, as the polypeptide progresses from a nearly random set of structures to the unique, compact native protein. An additional motivation for kinetic studies is their relation to the evolution of protein sequences. Evolution preserves protein sequences that correspond to structures with functions that are important to the organism. Theoretical studies by Wolynes and coworkers (2) have suggested how rapid folding to the native structure is yet another evolutionary pressure.
The experimental investigation of the kinetics and mechanism of protein folding has been aided by several recent theoretical and technological advances. The theoretical advances include analytical approaches (2–4), simulations of simplified representations of proteins (2, 5–8), and allatom molecular dynamics calculations (9–11). This work has painted a comprehensive picture of possible general mechanisms and has provided a framework for experimentalists to think more clearly about the problem. It also has helped define questions, design new experiments, and interpret experimental results. Important technological advances include the availability of a great variety of materials from protein engineering and peptide synthesis, the development of more rapid kinetic methods (12, 13), and increased computer power. The combination of these advances now permits the development of an “aufbau” approach to protein folding. This approach starts with the investigation of isolated secondary structural elements: αhelices, βstructures, and loops. The relative simplicity of these elements should permit their mechanism of formation to be described in much greater detail than is possible for proteins. Such studies include the development of statistical mechanical models which quantitatively reproduce equilibrium populations and kinetic progress curves. Once the kinetics and mechanism of the elements are understood, it should be possible to investigate structures of increasing size and complexity.
We have begun to study secondary structural elements by using nanosecondresolved kinetic methods and statistical–mechanical modeling (14). The thermodynamic and kinetic behavior of the αhelix has been studied for more than 40 years (15–20). Only recently, however, have kinetic measurements been made on helices of size and composition comparable with those found in proteins (21–23). Also, early theoretical studies (16, 17) were limited by the lack of computer power, preventing the detailed modeling of experimental kinetic data on helix formation that is now possible (13). The experimental and theoretical study of the kinetics of loops and βstructures is a new subject. Jones et al. (24) and Hagen et al. (25, 26) used a nanosecond photochemical triggering method to study loop formation in cytochrome c by determining the diffusion limited rate for an intramolecular ligandbinding reaction. We also recently reported a thermodynamic and kinetic study of a βhairpin formed by the 16 Cterminal residues of streptococcal protein G B1 (Fig. 1) (27). This peptide had been shown to adopt the βhairpin conformation by Blanco et al. (29) using NMR spectroscopy. Our βhairpin experiments consisted of measuring the thermal unfolding curve for the 16residue peptide between 273 K and 363 K and measuring the relaxation kinetics following 15degree nanosecond laser temperature jumps to final temperatures ranging from 288 K to 328 K (27). The three principal experimental results from this study were: (i) the βhairpin peptide exhibits twostate behavior in both its equilibrium and kinetics; (ii) the apparent activation energy for the folding rate calculated from the twostate analysis is negative; and (iii) the rate of βhairpin formation is much (>10fold) slower than that of the αhelices that have been studied up to now in short peptides.
To explain these results, we used a simple statistical mechanical model which was only briefly described (27). Here, we present a detailed description of the model, test one of its major approximations, and use the model to predict kinetic and equilibrium properties expected for other βhairpinforming peptides. We shall see that analysis of βhairpin thermodynamics and kinetics addresses many of the same issues that arise in considering the folding of a small protein.
Description of the Model.
Our objective has been to develop a model for protein secondary structure kinetics, which can be used to analyze experimental data and to predict new experiments. In this work, the model is applied to a βhairpin, but it also can be applied to helices and is readily adapted for more complex structures. We adopt a description, which uses pairs of φ,ψ dihedral angles to define the conformation of each molecule; the complete native structure is formed when all of the residues have native values for these angles. Formation of the native structure is opposed by the loss of conformational entropy and favored by the formation of stabilizing interactions, i.e., hydrogen bonds and hydrophobic interactions (Fig. 1). The model postulates that two groups interact only when all of the dihedral angles of the sequence connecting them are native. This restriction considerably simplifies the model by identifying threedimensional structures with sequences of peptide bond conformations.
A second simplifying step is to consider only two conformations for the backbone dihedral angles, native and nonnative (in a spirit similar to the “correct” and “incorrect” parameter of the Zwanzig model; ref. 30). The nonnative conformation of a dihedral angle pair is not a unique conformation but is the set of all conformations that are incompatible with the native structure. An additional feature of the model is that pairs of φ,ψ dihedral angles are assumed to rotate between native and nonnative values simultaneously.† We chose the dihedral angles ψ of residue i and φ of residue i + 1 (Fig. 2) so that the peptide bond, rather than the residue, is the conformational unit. Formation of a backbone–backbone hydrogen bond is therefore associated with the transformation of one pair of ψ_{i}, φ_{i+1} angles in each β strand from nonnative to native values.
In our thermodynamic description of the βhairpin, we consider only three factors. These are the stabilizing effect of the hydrogen bonds between the backbone carbonyl and amide of the N and Cterminal β strands, the stabilizing effect of the three hydrophobic interactions among the four side chains of the hydrophobic cluster (Fig. 1), and the destabilizing effect of the loss of conformational entropy when fixing pairs of dihedral angles in the native hairpin conformation. Nonnative interactions, such as wrong hydrogen bonds or hydrophobic interactions, are ignored. We also ignore electrostatic interactions among the charged side chains and chain termini (their importance could be assessed by experiments on the ionic strength dependence of the equilibrium and kinetics which have not yet been performed). Each thermodynamic factor is considered to be homogeneous, i.e., independent of side chain and position in the native structure. We assume that the free energies of formation for each of the three hydrophobic interactions, ΔG_{sc}, are identical. Each of the backbone–backbone hydrogen bonds, including the one in the turn region, is assumed to have the same free energy, ΔG_{hb}. The conformational entropy loss for the strand and turn regions also is assumed to be the same (ΔS_{conf}), which is equivalent to assuming that the residues in the turn have a propensity for this conformation equal to the propensity of the strand residues to be in a strand conformation. To further reduce the number of parameters, we assume that the hydrogen bond is purely enthalpic, i.e., ΔG_{hb} = ΔH_{hb} and that the hydrophobic interactions are temperatureindependent over the temperature range studied.
For the elementary kinetic steps (motion of individual dihedral angle pairs), we choose a transition state that can be described in terms of the equilibrium thermodynamic parameters. It is natural to assume that there is an entropy barrier to forming a native dihedral angle pair, so we equate the entropy of activation to the equilibrium entropy loss. For some steps, native dihedral angle pair formation is not associated with stabilizing interactions, whereas in others it is associated with the formation of hydrogen bonds or both hydrogen bonds and hydrophobic interactions. We assume that all native interactions are broken in the transition state. Also, we include the possibility that these steps have an activation barrier, E_{o}, in addition to the barriers imposed by the equilibrium free energy changes.
We must next decide how to treat the temperature dependence of the prefactor for these kinetic steps because it will have a significant effect on the height of the potential energy barrier required to fit the kinetic data. In investigating the viscositydependence of the conformational relaxation rate of myoglobin, Ansari et al. (31, 32) found that the data could be wellrepresented by a preexponential factor proportional to 1/(σ + η), where η is the solvent viscosity and σ is the contribution to the effective friction from interacting protein atoms (4 cP in myoglobin). A much greater fraction of the βhairpin peptide atoms interact with solvent, so we expect σ to be smaller. Simulations of βhairpin formation by Klimov and Thirumalai (33) suggest a 1/η dependence (σ = 0) so, in the absence of direct experimental data, we use a prefactor proportional to 1/η. The net result is that the model is completely defined by only five parameters—three equilibrium parameters, ΔH_{hb}, ΔG_{sc}, and ΔS_{conf}, and the two kinetic parameters, k_{o} (T_{o}), the preexponential factor at the reference temperature, and an activation energy, E_{o}.
A final, major simplifying feature in the model is the single sequence approximation first used by Schellman (34) in describing the helixcoil equilibrium and recently by us in describing helixcoil kinetics of a 21residue peptide (23). In the single sequence approximation, only species with a contiguous run of native peptide bonds are considered. All other structures are ignored. For the βhairpin peptide, which has 16 residues (15 peptide bonds), there are 2^{15} (= 32,768) possible molecular conformations. The single sequence approximation reduces this number to 121. In helixcoil theory, the justification for the single sequence approximation is the expectation that for short polypeptides there is a low probability of nucleating more than one stretch of helix in any individual molecule. For the βhairpin, we give the justification a posteriori by comparing with a more complete model in which the approximation is not made.
Partition Functions.
The nonnative conformation of the peptide bond (coil, c) is taken as the reference state and assigned a weight of 1. The weight of a peptide bond in the native conformation (hairpin, h) is exp(ΔS_{conf}/R), and the weight for a single stretch of j contiguous native peptide bonds, starting with peptide bond i [i.e., the ψ of residue i and φ of residue i + 1 (Fig. 2)], is: 1 where p is the number of backbone–backbone hydrogen bonds and q is the number of side–chainhydrophobic interactions in the native stretch. In this model, there are 2^{15} conformations for the 16residue hairpin arising from all of the possible combinations of hs and cs. The weight of each of these conformations is simply the product of the weights of each of the native stretches that it contains, and the partition function is the sum of the 2^{15} weights.
The model can be greatly simplified by considering only those species which contain a single stretch of native peptide bonds (the “standard” single sequence approximation). This simplification results in a model with only 121 species with the partition function: 2 where n + 1 is the total number of residues (16 in this βhairpin). The equilibrium probability of the allcoil conformation is P_{0,0} = 1/Q and the equilibrium probability for all other conformations is P_{j,i} = w_{j,i}/Q.
To test the accuracy of this standard single sequence approximation, we compared the equilibrium curves of the model with and without this approximation in Fig. 3. The approximation significantly overestimates the fraction of folded hairpin. This problem arises because the standard single sequence approximation does not properly account for the entropy of the system, as has been discussed by Qian and Schellman (35) for the helixcoil transition. The population of each of the 32,647 ignored species [such as cchcchcccchcccc with a weight of exp(3ΔS/R)] is quite small, but because their number is large, their contribution to the entropy of the system is significant. In particular, most of the ignored species do not contain significant hairpin structure, and ignoring them underestimates the stability of the unfolded hairpin. The number of species ignored by the standard single sequence approximation grows geometrically with peptide length, precluding its application to molecules of different length.
The underestimation of the entropy can be minimized by defining a “coil” state that includes not only the all c species (ccccccccccccccc) but also all the possible combinations of h and c peptide bonds that do not have just one single native stretch. For all those conformations, we ignore native interactions (even for a species such as ccchhhhhhhhhchc, which has the backbone conformation of the βturn as well as residues Y45 and F52 in position to make a hydrophobic interaction). The weight of the coil state now becomes: 3 where the second term eliminates the contribution to the coil state by conformations with a single stretch of native peptide bonds (see Eq. 1). The partition function in this “modified” single sequence approximation is: 4
Rate Equations.
To transform the equilibrium description of the model with the modified single sequence approximation into a kinetic model, we begin by assuming that conformations are connected if they can be interconverted by single h → c or c → h transitions. Species that contain a single stretch of native peptide bonds are connected to those other species that contain one more or one less native peptide bond at either end of the stretch. The rate constant for adding native peptide bond i to a native stretch, k_{i}^{+}, is given by: 5 where k_{o} is the preexponential factor at the reference temperature T_{o}, η_{o} is the solvent viscosity at T_{o}, and E_{o} is the activation energy for rotation of the peptide bond. The rate constants for removing native peptides bonds i or i + j − 1 from a native stretch of length j that starts at residue i (Fig. 2) are given by: 6 It is less straightforward to treat the contribution to the overall kinetics of the system of those additional species that now have been included in the coil state. For example, a coil conformation such as cchhhchhccccccc can convert to a single sequence conformation cchhhhhhccccccc by a single c → h transition. We assume that the rate for this process is equal to k_{6}^{+} (Eq. 5) times the probability of finding this particular conformation within the coil state (i.e., exp(5ΔS_{conf}/R)/w_{0,0} for the above example). We then can define an overall rate that is the summation of the rates for all possible transitions between the coil state and each conformation with a single native stretch. The overall rate for going from the coil state to a conformation with a stretch of j native peptide bonds starting at residue i is given by: 7 The overall reverse rate is given by: 8 where: Using these rates (Eqs. 5–8), the population of the 121molecular species of the model as a function of time is described by the following set of master equations: 9 where Despite its complexity, this treatment of the kinetics maintains detailed balance. Moreover, it implicitly includes all the kinetic connections involving single h → c or c → h transitions for each of the 120single sequence species without increasing the size of the rate matrix. The physical description in this approximation is, however, somewhat artificial. For example, our definition of the coil species requires that a c → h transition which does not occur at the end of a native stretch (such as ccchhhhhhhhcccc → ccchhhhhhhhchcc), transforms the molecule back to the coil state instead of closer to the fully formed hairpin. However, an additional single transition (ccchhhhhhhhchcc → ccchhhhhhhhhhcc) returns the molecule to a more complete hairpin conformation.
Test of Modified Single Sequence Approximation.
We tested the modified single sequence approximation by comparing it with a “complete” model that considers all 2^{n} (= 2^{15} = 32,678) possible conformations explicitly. To perform the test, we fit the experimental data with the modified single sequence approximation model to obtain parameters that were then used in simulations using the complete model. The fit and simulations of the equilibrium data are shown in Fig. 4a. The equilibrium description is rather similar for both models, in contrast to the standard single sequence approximation (Fig. 3). This result confirms our interpretation that underestimation of the entropy of the unfolded ensemble is the main deficiency in the standard single sequence approximation. In the modified single sequence approximation, however, there is a small overestimation of the fraction of unfolded hairpin. The major contribution to this difference is the small subset of species that has significant βhairpin structure (including stabilizing interactions) but are counted as species of the coil state (which have no stabilizing interactions) in the modified single sequence approximation.
We also tested the kinetic description with simulations carried out with the complete model, in which there are n2^{n} (= 491, 520) possible kinetic transitions. Only 450 of those are explicitly included in the modified single sequence approximation model. The fitting to the kinetic experiments with this model was performed by floating its five parameters to produce the best leastsquares fit to the observed progress curves for the nine experimental temperature jumps between 288 and 328 K (Fig. 4c). This was carried out using the equilibrium populations at the initial temperatures (before the Tjump) and integrating the rate equations (Eq. 9) by using rate constants evaluated at the final temperatures (after the Tjump). These parameters were then used in kinetic simulations with the complete model. The rate matrix for this model was constructed using an automatic patternmatching algorithm (E. R. Henry, unpublished data).‡
The results for the two models are very similar, with fluorescence progress curves that can be represented as a biexponential process in each case. There is initially a small amplitude phase, corresponding to very rapid reequilibration among conformations in the globalfree energy minima of the folded state, followed by a slower largeamplitude phase, corresponding to crossing of the free energy barrier separating the folded and unfolded states (see Fig. 5 and below). Overall, the agreement between the two models must be considered very good and justifies the use of the modified single sequence approximation. There are, however, significant differences, and the relaxation rates for the major phase (the only one detected experimentally) are about a factor of three faster in the complete model (Figs. 4 b and c). This effect is produced because the modified single sequence approximation ignores the stabilizing interactions in the rates connecting conformations included in the coil state with the conformations in the folded state (with a stretch of seven or more native peptide bonds). For example, the transition cccchhhhhhchccc → cccchhhhhhhhccc is less probable in the simpler model because ignoring the two hydrogen bonds of the starting conformation lowers its population by a factor of 25 [= exp(−2ΔH_{hb}/RT)].
Predictions for Other βHairpins.
An important consequence of having a statistical mechanical model for βhairpin formation is that it can be used to make specific predictions that can be tested experimentally. A useful way of examining the results of the model is to consider the free energy as a function of the fraction of native peptide bonds, its natural reaction coordinate (Figs. 5a and 6a). The model postulates that formation of a βhairpin in the absence of side–chain interactions is continuously uphill in free energy because backbone hydrogen bonds do not compensate for the loss in conformational entropy of forming native peptide bonds in both βstrands. Side–chain interactions are, therefore, necessary for the stability of the hairpin and determine the position and height of the free energy barrier for hairpin formation. This hairpin is stabilized by a cluster of four hydrophobic side chains (W43, Y45, F52, and V54), making three hydrophobic interactions (Fig. 1). Based on our model, the folding free energy barrier for this hairpin is crossed when the seven central peptide bonds become native and the first hydrophobic interaction (between residues Y45 and F52) is formed.
Many of the predictions of the model are immediately apparent upon examination of the free energy profile. The existence of two global minima separated by a significant free energy barrier (Fig. 5a) explains the twostate behavior and exponential kinetics. The species at the barrier maximum has two backbone–backbone hydrogen bonds and therefore a lower energy than the coil state, explaining the apparent negative activation energy in the twostate analysis (assuming a simple Arrhenius expression for the rate constants with a temperatureindependent prefactor). The global minimum on the folded side of the free energy barrier consists of several molecular conformations, with the species at the lowest free energy having the intact hydrophobic cluster but not the maximum number of backbone–backbone hydrogen bonds and native peptide bonds. This result could explain why the population of the hydrophobic cluster obtained by fitting fluorescence data is higher than the fraction of native dihedral angles estimated by NMR (29).
An interesting prediction of this model of βhairpin formation is that local and longrange interactions have very different effects on the free energy surface of the hairpin and, therefore, on its equilibrium and kinetic properties. To illustrate this point, we have performed simulations with two variants. In one of these variants, we include a side–chain interaction, which could result from a favorable electrostatic interaction between D47 and K50, which stabilizes the βturn by 1 kcal/mol. In the other variant, a similar interaction is introduced between the first and last residues in the hairpin by mutating “in machina ” glycine 41 to arginine. This computational experiment is similar in spirit to the protein engineering approach to folding kinetics (41–43). When positioned in the βturn, the interaction is local (between residues i and i + 3), and it significantly affects both the thermodynamics and kinetics of hairpin formation by lowering the free energy of all states, which contain native interactions (Fig. 5a). Both the population of species with the hydrophobic cluster and the fraction of hydrogen bonds increase at all temperatures, and the T_{m} increases by ∼20 K (Fig. 5b). The folding rate is accelerated by about a factor of four, whereas the unfolding rate is slightly decelerated (Fig. 5c). Because the peak of the free energy barrier stays at the same position along the reaction coordinate, the change in rates results simply from the change in the barrier height, as is commonly assumed in interpreting the effects of single residue perturbations in protein folding. When the interaction is introduced between the end residues, it is long range (i,i + 15) and its effects on the folding properties are rather insignificant. The T_{m} changes by only ∼2 K, and the change in rates is very small, with the largest change in the unfolding rate. Thus, the simulations suggest that, in hairpins, the interactions closest to the βturn exert the largest effect on the folding rate; interactions between the ends of the strands may stabilize the hairpin structure but have very little effect on the folding rate.
Another important point raised by our model of βhairpin formation is that the shape of the free energy barrier and the position of its maximum along the reaction coordinate are determined by a delicate balance between the loss in conformational entropy and stabilization from side–chain interactions. To address this point, we have simulated two variants of the original βhairpin: a hairpin with the hydrophobic cluster placed one residue closer to the center of the molecule (W44, Y46, F51, V53) and another one with the hydrophobic cluster one residue closer to the ends (W42, Y44, F53, V55). The effects of these changes on the equilibrium and kinetic properties are shown in Figs. 6 b and c, respectively. Moving the hydrophobic cluster one residue in either direction does not modify the interaction energies of the hairpin; however, the model predicts a dramatic effect on its free energy profile (Fig. 6a). If the cluster is moved closer to the βturn, both the minimum in the folded ensemble and the top of the free energy barrier are shifted toward less structure (closer to the unfolded ensemble). This is accompanied by an increase in stability, an acceleration of the folding rate, and almost no change in the unfolding rate. If the cluster is moved one residue in the opposite direction (toward the ends), the stability is decreased and the folding rate decreases, and there is only a small change in the unfolding rate. Displacements of the top of the free energy barrier have been reported in folding experiments on small proteins (44). The model indicates that, for a βhairpin, the top of the free energy barrier is simply determined by the position of the stabilizing side chains in the sequence.
Caveats.
One criticism of the model that we have presented is that only native interactions are considered. This excludes the possibility, for example, of forming a turn at additional positions in the sequence, which would result in nonnative hydrogen bonds and nonnative hydrophobic interactions. There is no evidence in the NMR data of significant population of other hairpin conformations (29). Nonnative interactions also can affect the kinetics in the same two ways as in proteins (2). They can produce local minima in the energy landscape, which can result in the population of intermediate structures at equilibrium, or they can produce transient trapping of misfolded structures, which are not present at equilibrium. We have not yet found any evidence for equilibrium intermediates in the folding of the βhairpin. For transient trapping to be observable as a separate kinetic phase, the residence time in the trapped state must be longer than the relaxation time for the overall hairpincoil transition. Transient trapping does not appear to be occurring in this βhairpin because the progress curves at all temperatures can be wellfit with a single exponential function (27).
Another criticism of the model is that there are no native backbone or side–chain interactions between two residues unless the peptide bonds of all intervening residues have the native conformation. This postulate excludes the possibility, for example, of initiation by forming the hydrophobic cluster followed by zipping up of the hydrogen bonds. The transition state in this mechanism would be a ∼10residue loop. One is tempted by this mechanism because of the close correspondence of the βhairpin relaxation time and the time of ∼1 μs estimated by Hagen et al. (25, 26) to form a 10residue loop. A 10residue loop is predicted by Thirumalai (45) to be the most probable loop size in proteins, longer loops being less probable because of the larger entropy loss and shorter loops because of chain stiffness. One could possibly distinguish between the two mechanisms by measuring the kinetics for a βhairpin in which the hydrophobic cluster is moved closer to the βturn. Our model predicts that the rate of formation should speed up because the transition state now occurs earlier along the reaction coordinate (Fig. 6a), whereas the loop model would predict a slower rate of formation. This consideration was in fact one of the motivations for the predictions discussed above.
The most convincing test of the model will of course come from measurements on other βhairpin peptides. Another approach to both testing and refining the model is to examine the results of simulations. Allatom molecular dynamics simulations of temperature jump kinetic experiments (at experimental temperatures) may be feasible in the near future (10, 11). In the meantime, it should be useful to examine the results of Langevin simulations of simplified representations of the peptide (33). Because large numbers of sufficiently long trajectories are possible with this method, kinetic progress curves actually can be simulated. Examination of these trajectories might reveal dominant mechanisms and structural species that must be included in a kinetic model. The results will, however, depend critically on the choice of potential functions.
Is the model unnecessarily complex? Could we use a much simpler model, even a twostate model? The main problem with a twostate model is that it has little predictive value. In a twostate model, one would postulate a transition state structure or, as in the case of proteins, try to determine the transition state by structural perturbation experiments. The examples in Fig. 6 show that small structural perturbations lead to changes in the transition state and would therefore also result in incorrect predictions for the change in rates. Nevertheless, the model could be simplified. If we consider only residue–residue interactions, then in the single sequence approximation, the model reduces to an eightstate model. In such a model, quartets of dihedral angles change simultaneously in single kinetics steps. This model can explain the experimental data, as well as predict the properties of other βhairpins. It is, however, not straightforward to extend a model based on interactions to more complex structures because of the difficulty in defining rules that specify the rates of all elementary kinetic steps in terms of just a few parameters.
Acknowledgments
We thank Attila Szabo and Peter Wolynes for helpful comments on the manuscript.
Footnotes

↵* To whom reprint requests may be addressed: email: vmunoz{at}helix.nih.gov or eaton{at}helix.nih.gov.

This paper was presented at the colloquium “Computational Biomolecular Science,” organized by Russell Doolittle, J. Andrew McCammon, and Peter G. Wolynes, held September 11–13, 1997, sponsored by the National Academy of Sciences at the Arnold and Mabel Beckman Center in Irvine, CA.

↵† When pairs of dihedral angles are used instead of single dihedral angles, the specification of a pair of angles produces a problem in phasing between the loss of entropy and the compensating decrease in interaction free energy. Either choice of φ,ψ pairs represents a compromise. This can be illustrated by considering the formation of a sixresidue βhairpin with a side–chain interaction between residues two and five. To form the backbone–backbone hydrogen bond requires native values for four dihedral angles, φ_{3},ψ_{3},φ_{4},ψ_{4}. If we were only concerned with hydrogen bond formation, as in helixcoil theory for homopolypeptides, then the natural choice for the dihedral angle pairs would be the φ and ψ associated with the same residue—in this case the two pairs φ_{3},ψ_{3} and φ_{4},ψ_{4}. With this choice, however, formation of the two to fiveside–chain interaction requires that eight dihedral angles assume native values—when only six, i.e., ψ_{2},φ_{3},ψ_{3},φ_{4},ψ_{4}, and φ_{5}, actually are required. So, in choosing ψ_{i},φ_{i+1} instead of φ_{i},ψ_{i} pairs, we overestimate the loss in entropy associated with formation of the first hydrogen bond, in favor of accurately representing the compensation between entropy loss and formation of side–chain interactions in subsequent steps.

↵‡ The system of equations is stiff and was integrated using an iterative multistep backward differentiation formula method (37), as implemented in the cvode package (36, 38). This algorithm requires the solution of a set of nonlinear algebraic equations by Newton iteration at each time step. Each Newton iteration in turn requires solving an NxN linear system AΔP = residual, where the matrix A is derived from the rate matrix K. For n = 32,768, this problem is rather too large to solve using standard methods (39). However, the matrix A is sparse, containing only ≈500,000 nonzero elements of a possible 10^{9}. Therefore, an iterative generalized minimal residual method (40) appropriate for large sparse linear systems, as implemented in the cvode package (36, 38), was used. The performance of the algorithm was improved dramatically in this application by Jacobi (diagonal) preconditioning or very simple blockdiagonal preconditioning (40).
References
 ↵
 Anfinsen C B
 ↵

 Bryngelson J D,
 Wolynes P G
 ↵
 Elber R
 Orland H,
 Garel T,
 Thirumalai D T
 ↵
 ↵
 ↵
 Guo Z,
 Brooks C B,
 Boczko E M
 ↵
 ↵
 Lazaridis T,
 Karplus M
 ↵
 ↵
Gray, H. B. & Valentine, J. S. (1998) Acc. Chem. Res., in press.
 ↵
Eaton, W. A., Hofrichter, J., Muñoz, V. & Thompson, P. A. (1998) Accts. Chem. Res., in press.
 ↵
 Zimm B,
 Doty P,
 Iso K
 ↵
 ↵
 Poland D,
 Scheraga H A
 ↵
 ↵
 ↵
 ↵
 Jones C M,
 Henry E R,
 Hu Y,
 Chan CK,
 Luck S D,
 Bhuyan A K,
 Roder H,
 Hofrichter J,
 Eaton W A
 ↵
 Hagen S J,
 Hofrichter J,
 Szabo A,
 Eaton W A
 ↵
 Hagen S J,
 Hofrichter J,
 Eaton W A
 ↵
 ↵
 Gronenborn A M,
 Filpula D R,
 Essig N Z,
 Achari A,
 Whitlow M,
 Wingfield P T,
 Clore G M
 ↵
 ↵
 Zwanzig R
 ↵
 Ansari A,
 Jones C M,
 Henry E R,
 Hofrichter J,
 Eaton W A
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Hairer E,
 Wanner G
 ↵
 Cohen S C,
 Hindmarsh A C
 ↵
 Lawson C L,
 Hanson R J
 ↵
 Barrett R,
 Berry M,
 Chan T F,
 Demmel J,
 Donato J M,
 Dongarra J,
 Eijkhout V,
 Pozo R,
 Romine C,
 Van der Vorst H
 ↵
 ↵
 ↵
 ↵
 Camacho C J,
 Thirumalai D