- Email: [email protected]
Progressing from folding trajectories to transition state ensemble in proteins
Chemical Physics 307 (2004) 251–258 www.elsevier.com/locate/chemphys
Progressing from folding trajectories to transition state ensemble in proteins D.K. Klimov a
a,*
, D. Thirumalai
b
Bioinformatics and Computational Biology Program, School of Computational Sciences, George Mason University, Manassas, VA 20110, USA b Department of Chemistry and Biochemistry, Institute for Physical Science and Technology, University of Maryland, Bldg 085, College Park, MD 20742, USA Received 5 March 2004; accepted 25 June 2004 Available online 27 August 2004
Abstract A general procedure for determining the nature of the transition state ensemble (TSE) in proteins is used to probe the heterogeneity of TSE in two models of a/b proteins. The sequences, which differ only in the chain connectivity, are modeled using a-carbon representation. The structures in the TSE are determined using a progress variable clustering (PVC) algorithm that does not depend on the choice of the reaction coordinate. Folding simulations starting from these transition state structures partition with nearly equal probability into the native state or unfolded state. Our fully automated PVC algorithm, which utilizes only the time series in the folding trajectories, is also computationally efficient. The results for the model a/b proteins show that the TSE can be divided into a small number of clusters (two or three) consisting of structurally related conformations. For both sequences the conformations in the dominant transition cluster are compact and native-like with the secondary structure elements fully formed. Although none of tertiary contacts is completely formed in the transition region, few occur with high probability. This suggests that there are multiple (more than one) nuclei that connect the unfolded and native states. Overstabilizing the high probability nucleation contacts retards folding rates. The results for both the models also suggest that, if the transition region is close to the native state, then permutation of the secondary structural elements does not alter the nature of the transition state. Ó 2004 Elsevier B.V. All rights reserved. Keywords: Protein folding; Transition state; Multiple folding nuclei; Progress variable clustering method
1. Introduction The transition from unfolded conformations to the folded state of proteins and RNA involves multiple pathways, in which a number of non-covalent interactions form and break at early times. For two state proteins, near transition state (TS) region a fraction of relatively stable (‘‘nucleating’’) contacts form. Upon formation of the nucleating contacts, which is a mixture of both long and short range contacts (range here refers to the distance along the sequence), the transition to the *
Corresponding author. Tel.: +3014054786; fax: +3013149404. E-mail address: [email protected] (D.K. Klimov).
0301-0104/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.chemphys.2004.06.071
native state occurs with overwhelming probability [1]. The consolidation of the entire structure immediately after the nucleating event is accompanied by the global chain compaction, so that folding and collapse are nearly synchronous. Because of the complexity of protein folding [2–8] it is difficult to obtain a simple reaction coordinate to describe the kinetics of protein folding. Given the multidimensional nature of the free energy landscape physical arguments are often used to obtain a low dimensional pseudo-reaction coordinates [9–11]. For example, for simplified protein models a fraction of native contacts, Q, which measures the assembly of native interactions, has been suggested as a suitable reaction coordinate [9,12]. Several studies demonstrated
252
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
the utility of Q in estimating the free energy barriers for folding as well as the approximate location of the transition state for small two-state folding proteins [9]. In general, it has been argued that Q is an appropriate reaction coordinate provided the protein sequence is highly optimized to fold to the native structure. Unlike in simple chemical reactions involving bond breakage or simple conformational transitions (boat to chair transition in cyclohexane, for example), the very definition of the TS (more precisely, the transition region) is ambiguous. Theoretical arguments suggest that, typically, the transition region is broad with superimposed structure so that it is fruitful to visualize a transition state ensemble (TSE) rather a unique TS for protein folding [13,14]. Such a plasticity in the TSE also rationalizes the ability of many proteins to tolerate large number of mutations without compromising the native fold. Experimental determination of the TS is based on /-value analysis pioneered by Fersht [15,16]. Although the interpretation of /-values is not always straightforward, those close unity indicate that a particular residue is highly structured in the TS. Combined experimental and simulation efforts have given considerable insight into the nature of TSE in a number of proteins. Nevertheless, the straightforward implementation of /-value analysis has been of limited use in the simulations of protein models including the all-atom ones due to considerable computational efforts in reliable determination of thermodynamic and kinetic properties of wild-type proteins [11]. Moreover, the determination of free energies requires the use of suitable order parameter, which can be ambiguous [17]. The computational challenge is to establish the nature of the folding TSE using only the time series generated by suitable dynamics. In other words, given the dynamics of all the particles starting from initially denatured state till the first passage time is reached for a large number of events how can one determine the TSE without making any assumption about the underlying reaction coordinate? Because all-atom molecular dynamics simulations cannot be used to generate a single trajectory for any protein till it reaches the native state starting from unfolded conformation, most of the analysis of TSE has been done using coarse-grained models. In our previous studies we developed a method based on clustering the structures in trajectories and computing TSE by using the ‘‘progress variable’’ d = t/s1i, where s1i is the first passage time for the trajectory i and t s1i [10,14]. Because this method relies both on the progress variable d and clustering of structures, we refer to it as PVC method. In implementing PVC method, we defined the TS structure for a given trajectory as the one which contains the minimal set of stable native contacts, the formation of which ensures rapid passage to the native structure. In practice, we seek a collection of conformations, a putative TSE, which corresponds to sharp rise in
the fraction of native contacts as a function of d. By computing these structures for a large number of independent folding transitions, a meaningful description of the folding TSE may be obtained. An alternative definition of TSE follows from the work of Klosek et al. [18]. In this approach the TS region is taken to be the locus of structures (stochastic separatrix) that have equal probability to lead either to unfolded U or native N states, i.e., the fraction of folding (‘‘forward’’) trajectories pfold 0.5. This idea was first numerically implemented by Du et al. [19] to describe folding kinetics in lattice models without side chains. It has been subsequently applied to more complex lattice models with side chains and other systems [10,20]. Although neither the PVC method nor the pfold technique relies on the choice of a reaction coordinate, they are both subject to certain limitations. The first provides approximate characteristics of the TSE, because its position along the folding pathway is determined from the average folding trajectory. The PVC method is computationally efficient, because it utilizes only the folding trajectories without the need for additional simulations to compute TSE. The stochastic separatrix (pfold) method assumes that at T TF a protein, upon crossing TS, subsequently folds to the native state or unfolds with approximately equal probability [19]. Therefore, if one initiates folding simulations with TSE conformations, the fraction of molecules pfold that finds the native state prior to unfolding should be 0.5. Obvious advantage of pfold method is that the saddle point on the free energy landscape is located by using a physically reasonable (but by no means unique) definition of TS. Yet the obvious computational inefficiency of the pfold method severely limits its practical implementation. In addition, there is some ambiguity in characterizing the unfolded basin of attraction that can compromise the numerical implementation. The purpose of this paper is to compare the TSEs determined by both methods for the coarse-grained models with a/b protein motif. Both coarse-grained Ca-models of a/b proteins differ only in the connectivity of the secondary structural elements (SSEs). For both of those we find that the two methods that use completely different definitions of the TSE give a consistent picture. Although there are appreciable differences in the heterogeneity of folding pathways in the two models the nature of the TSE are quite similar. This surprising result is rationalized in terms of the nearness of the TSE to the native state. 2. Methods 2.1. A coarse-grained model and simulation details for a/b-protein We used a coarse-grained off-lattice Ca-carbon model for an a/b protein. The model sequence consists of
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
N = 38 connected beads of three types, hydrophobic B, hydrophilic L, and neutral N. The potential energy includes contributions from bond-length and bond-angle potentials, the dihedral angle potential, and non-bonded potential described in detail elsewhere [21,22]. The average separation between Ca-carbons along the sequence is ˚ . A sequence is determined by the connectivity a = 3.8 A of residues and the set of random factors k in the nonbonded potential [21]. Two sequences, A1 and A2, were designed to adopt a/b-motifs in their native state, which contains two b-strands (s1 and s2) and a-helix (h) connected by turns (Fig. 1). The a/b topology is stabilized by proper dihedral angle potentials and distribution of hydrophobic residues and k factors achieved using Monte-Carlo sequence design algorithm [23,24]. We used Langevin dynamics (LD) to simulate the folding of the a/b protein. The thermodynamics of folding is computed in the underdamped limit by applying the velocity form of Verlet algorithm [21]. Multiple histogram method was used to compute several probes, such as the fraction of native contacts Q, the fractions of native contacts between SSEs Qij (i,j = s1, s2, h), and the radius of gyration Rg as a function of temperature T. The folding kinetics is obtained from LD simulations in the overdamped limit (approximately at water viscosity) at the temperature Ts . TF, where TF is the folding temperature. The folding transition temperature is determined using standard methods [10]. The mean first passage time sfp is computed by averaging over 400 trajectories. 2.2. Analysis of folding TSE The structures of the TSE are determined using the clustering (pattern recognition) algorithm [10,13,14]. The conformations from folding trajectories generated at Ts are recorded over regular intervals ds = 120 ps(1 + [s1i/600 ns]), where the first passage time s1i expressed in ns corresponds to the first instance of assembly of all native contacts. We tested different values of ds
253
to ensure that qualitatively similar results are obtained. In total, we used 100 folding trajectories. The recorded conformation k is represented by a vector Ck with the components given by the distances between the residues forming native contacts rmn. The vectors Ck (k = 1,. . .,Ns,i, where Ns,i is the number of conformational snapshots for the trajectory i) are clustered to filter out fast structural oscillations [14]. Clustering procedure for a trajectory i requires choosing the values of two parameters, the cluster radius Rc,i and the tolerance for native contacts Dc,i, which defines the maximum deviation of rmn from the native value r0mn . Four criteria are used to select Rc,i and Dc,i: (i) the ‘‘native’’ cluster, which includes the native conformation, should have all native contacts formed; (ii) there shall be no other clusters beside native, in which all native contacts are formed; (iii) the number of structures in the native cluster must be small ([3–4); (iv) the native cluster must first appear close to s1i (e.g., d J 0.99). The optimal values of Rc,i and Dc,i are localized in the upper right corner of the plot in Fig. 2. This choice for Rc,i and Dc,i ensures that the number of clusters Nc,i satisfies the condition 1 Nc,i Ns,i, making the clustering procedure meaningful. The stability of contacts until s1i is determined with a certain tolerance, which permits short-lived disruptions. Following our earlier studies [14,25], we assume that the ‘‘nucleation contacts’’ at a given d constitute the minimal set of native contacts that remain stable until s1i. To establish the location of the putative TSE region, we calculated the kinetic probability Pq(d) averaged over all folding trajectories that a contact q is part of the minimal set of stable native contacts at a given d. From this, we obtain the average (over all native contacts) probability P(d) for a native contact to be part of a minimal set of stable contacts. The value of dTS, at which dP/ dd begins steep rise, is taken to be the TS region. The entire procedure of locating TSE, starting from the clustering of folding trajectories to computing P(d), does not depend on the protein model and is fully automated.
Fig. 1. Ribbon diagrams of native conformations for sequences A1 and A2 visualized using MolMol [35]. Although both native conformations for A1 and A2 have an a/b motif, their topologies are different. A1 adopts strand–strand–helix (b–b–a) fold, whereas the native state of A2 is strand– helix–strand (b–a–b).
254
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
Fig. 2. Diagram of clustering parameters values (Rc,i,Dc,i) meeting the four requirements for cluster procedure described in the text (shaded area). The optimal values of Rc,i and Dc,i correspond to the upper right corner of the shaded area. Efficient elimination of fast structural fluctuation is achieved when Rc,i is as large as possible. The data are plotted for a generic trajectory for sequence A2.
Thus, our method uses only the time series for a number of folding trajectories to determine TSE.
3. Results 3.1. Folding transition state ensemble The TS region has been located using the kinetic probability P(d) (see Section 2). The derivative dP/dd (and P(d), data not shown) remains close to zero for the most part of folding trajectories, implying the absence of stable native interactions (inset to Fig. 3(a)). The range of d, corresponding to the beginning of dramatic rise in dP/dd, is associated with the TS region. After tests of various d we concluded that TS in both sequences, A1 and A2, is located at dTS = 0.9 (inset to Fig. 3(a)). The particular value of dTS is not important as long as dTS is kept within a certain interval (0.85 [ dTS [ 0.95). However, at d < 0.85 a large portion of TS conformation is completely unstructured, whereas at larger d > 0.95 the entire TSE is compact and homogeneous. Propagation of stable native contacts Pq(d) for both sequences (Fig. 4) reveals that, in contrast to rapid assembly of local native helical contacts, stable interactions between SSEs emerge only relatively late in folding process (large d) following the formation of compact transient intermediates. Furthermore, stable tertiary contacts across the sequences form at roughly similar d. As a result the TS is compact and the folding nuclei are dispersed across the sequence [6,26]. These findings imply that the formation of stable native folded structure requires tertiary interactions [27]. Similar propaga-
tion of stable native interactions in both sequences implies that the differences in topological connectivity are not evident in the TSE. This results is a consequence of convergence of folding pathways prior to crossing the TSE. To get a better insight into the formation of native interactions in TSE we consider the probability Pn that the residue n is structured in the TSE (i.e., is engaged in nucleation contacts). From the profiles of Pn at several d, including dTS = 0.90 (Fig. 3) the following conclusions can be drawn: (a) The majority of residues in A1 and A2 are involved to a varying degree in the folding nuclei (with the exception of a few, mostly in turn regions, which do not make native interactions), but none of the contacts is fully structured (Pn = 1) [14,28,29]. As a result the nucleus is delocalized and the TS structures are compact [6,26]. (b) Although most of residues are involved in the folding nuclei to some extent, there are regions, which are more structured than the others. Due to topological connectivity, residues involved in nucleation are often located near the turns, connecting the SSE as predicted earlier [13]. This is seen for the residues in the A1 hairpin adjacent to the s1–s2 turn. In A2, the residues in the strand s1 and helix are more structured than those in s2, which suggests that the TS in A2 is more ‘‘polarized’’ than in A1. The interactions formed in the TS are expected to strongly affect folding times. To verify this expectation we selected in A1 eight TS hydrophobic contacts with high Pq(dTS). Four mutants were designed by increasing
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
(a)
(b) Fig. 3. Kinetic probabilities Pn for the residue n at three values of d close to the TS value dTS = 0.90, for A1 (a) and A2 (b). Pn, which estimates the degree to which a particular residue n is structured in the TS, shows that folding nuclei are diffused as the majority of residues participate to a varying degree in the TSE. The inset in (a) shows dP/ dd for A1 (see Section 2). The value of dTS = 0.90 marked by dashed vertical line corresponds to the beginning of dP/dd growth, which is associated with crossing TS. dP/dd for A2 shows similar behavior. Data are calculated using 100 folding trajectories.
the strength of these interactions by 0.63, 1.25, 1.88 and 2.50 kcal/mol (variation by about 8–30%). Compared to the WT the folding times of the first two mutants decrease (to 0.93sF and 0.89sF). However, further stabilization of TS interactions (last two mutants) impedes folding and the folding time increases to 1.05sF and 1.27sF. Therefore, the strength of TS interactions must be optimal as too strong TS contacts may lead to entanglement traps, which in turn slow down folding [14]. 3.2. Structural diversity of the TSE To assess structural diversity of TSE, we applied the clustering analysis to TSE structures [25]. The A1 TSE
255
can be divided into three distinct clusters. The first, TS1, is compact and native-like. The ratio of the radius of gyration of TS1 conformations to its native value Rg =RN g is 1.0. TS1 contains, on an average, about twothirds of native contacts (Q = 0.68). The secondary structures in TS1 are completely formed. Most of the tertiary interactions between SSE are formed with the hairpin s1–s2 being the most structured region. TS1 appears in 70% of all folding trajectories. The structures in the TS2 cluster are less compact ðRg =RN g 1:3Þ and structured than TS1. In TS2 the hairpin is mostly formed, while there are no significant interactions between the helix and strands. The secondary structure in TS2 is not as fully formed as in TS1 (the b-strand and a-helix contents are between 0.5 and 0.6). About 25% of folding trajectories pass through TS2. The structure of rare TS3 is similar to TS2, but differs in the location of the structured region. In TS3, the tertiary interactions are formed only between the strand s2 and h. There are two dominant TS clusters in A2, the first of which, TS1, is very similar to TS1 of A1. This TS is compact ðRg =RN g ¼ 1:0Þ and native-like (Q = 0.60). About 80% of unfolded conformations reach the native state through TS1, while the remaining 20% of molecules fold via the less structured TS2 ðQ ¼ 0:22 and Rg =RN g ¼ 1:4Þ. More importantly, the tertiary interactions between SSEs in TS2 are found only between the strand s1 and helix. 3.3. Application of pfold to TSE For two-state folders the TSE may be defined as a set of structures, from which the probabilities to reach native N or unfolded U states are approximately equal. Folding trajectories initiated with TSE conformations should partition into those, which reach N prior to unfolding (‘‘forward’’ trajectories) and those, which initially unfold. The fraction of ‘‘forward’’ trajectories, pfold, should be approximately 0.5. This requirement is physically reasonable operational definition of the TSE [19]. To assess whether the TSE identified by the PVC method is consistent with the pfold definition, we calculated pfold for the TS structures obtained at dTS = 0.90. In all, we recorded 129 putative TS conformations for A1 and 152 for A2 and selected every fifth conformation. Each such conformation serves as a starting point for M = 100 folding trajectories generated at T = Ts. Each trajectory either reaches the native state (Q = 1.0) or leads to U (i.e., when the number of tertiary contacts Cu < 5). The value of Cu was varied to ensure that the resulting pfold values do not change significantly. For the kth putative TS conformation we calculated the coefficient pkfold ¼ M kF =M, where M kF is the number of ‘‘forward’’ trajectories. Then pfold is taken to be an
256
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
(a)
(b) Fig. 4. The kinetic probabilities Pq as a function of d averaged over 100 trajectories for A1 (a) and A2 (b). Pq gives the probability that a native contact q belongs to a minimal set of stable native contacts at a given d which survives until the first passage time. The color codes to Pq are given on the right. The plots show a remarkable heterogeneity in the formation of native contacts. Helical q become stable very early in the folding process, while all other native contacts irrespective of their location stabilize close to reaching the native state, i.e., the TS occurs late in the folding process. Surprisingly, plots for A1 and A2 are qualitatively similar. The native contact indexes q are arranged according to the type of secondary structure and sequence separation.
average over all pkfold . The values of pfold for A1 and A2 are equal to 0.51 and 0.59, respectively. Thus, the descriptions of the TSE obtained by the stochastic separatrix method and the PVC method are consistent.
3.4. Chain connectivity and the nature of TSE The dominant transition state (TS1) for both the models have similar characteristics (Fig. 1) despite the fact
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
that A1 and A2 have different connectivity between the SSEs. Experimentally the variations in the chain connectivities can be achieved by circular permutation. Whether or not the TSE shifts upon a circular permutation seems to depend on the protein topology and the precise location of permutation. It has been suggested that the nature of the TSE does not change upon altering chain connectivity for CI2 [30,31], but does so dramatically for the SH3-domain of a-spectrin [32,33]. Our simulations suggest that, if the structures in the dominant TS resemble closely that of the native state, then circular permutation would have little effect on the TSE. In both the model sequences structures in the dominant TSs are nearly nativelike with fully formed secondary structures. Such a highly native-like structured TS is only possible, if earlier events lead to specific collapse. The similarity in the structural details of the dominant TSs of two models suggests that the TSE is crossed in the last stages of folding. Experimentally this would imply that if the Tanford b values are close to unity, as is the case in a class of cold-shock proteins [34], then the TSE would not alter appreciably upon permutation provided that circular permutation does not alter the fold itself. For this class of proteins we predict that the TSE would be relatively homogeneous for all the mutants.
4. Conclusions In this study we have used two independent methods for analyzing folding TSE. The stochastic separatrix (pfold) method, which provides a physically motivated description of TSE, confirms the location of TS established using PVC method. Through the analysis of hundreds of folding trajectories we show that dominant folding TS exists for both sequences, A1 and A2. Most molecules (70–80%) fold through this TS, which is compact and topologically resembles the native fold. In the dominant TS the SSEs are fully formed. This conclusion is consistent with the calculated time scales for collapse sc and folding sF. The ratio sF/sc is 2.7 for A1 and 2.4 for A2. These results together with the TS analysis show that the dominant folding TS (corresponding to the ratelimiting step in folding) is native-like, compact, and involves interactions (to a varying degree) between all SSEs. This picture of TSE resembles the notion of delocalized folding nuclei diffused across the native structure [6,9,26]. Because the rate-limiting step occurs late in the compact phase, the dominant folding TSs for A1 and A2 are quite similar, even though their native states are topologically different. Despite the existence of a dominant single TS, there is still certain diversity in the TSE as 20– 30% of molecules fold through alternative TSs. The structure of TS2 is polarized, because native interactions are predominantly formed near the N-terminal (the bhairpin in A1 or s1–h construct in A2).
257
This study shows that a combination of two methods based on d progress variable (a PVC method) and the notion of stochastic separatrix (pfold method) offers a powerful tool for locating and analyzing folding TSE. Because both the methods are general and do not rely on a choice of a folding reaction coordinate, they may be applied to variety of protein models, from simple lattice to all-atom representations. The computationally more efficient PVC method is easier to use.
Acknowledgment This work was supported in part by the National Science Foundation Grant No. NSF CHE-0209340.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
[16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]
P.G. Wolynes, Proc. Natl. Acad. Sci. USA 93 (1996) 14249. A.R. Fersht, Proc. Natl. Acad. Sci. USA 92 (1995) 10869. Z. Guo, D. Thirumalai, Biopolymers 36 (1995) 83. E. Shakhnovich, V. Abkevich, O. Ptitsyn, Nature 379 (1996) 96. A.R. Fersht, Curr. Opin. Struct. Biol. 7 (1997) 3. P.G. Wolynes, Proc. Natl. Acad. Sci. USA 94 (1997) 6170. J.N. Onuchic, Z.A. Luthey-Schulten, P.G. Wolynes, Ann. Rev. Phys. Chem. 48 (1997) 545. L. Mirny, E. Shakhnovich, Ann. Rev. Biophys. Biomol. Struct. 30 (2001) 361. N.D. Socci, J.N. Onuchic, P.G. Wolynes, J. Chem. Phys. 104 (1996) 5860. D.K. Klimov, D. Thirumalai, Prot. Struct. Funct. Genet. 43 (2001) 465. C. Clementi, H. Nymeyer, J.N. Onuchic, J. Mol. Biol. 298 (2000) 937. C. Clementi, P. Jennings, J. Onuchic, Proc. Natl. Acad. Sci. USA 97 (2000) 5871. Z. Guo, D. Thirumalai, Folding Des. 2 (1997) 377. D.K. Klimov, D. Thirumalai, J. Mol. Biol. 282 (1998) 471. A.R. Fersht, Structure and Mechanism in Protein Science: A Guide to Enzyme Catalysis and Protein Folding, Freeman, New York, 1998. A.R. Fersht, Curr. Opin. Struct. Biol. 5 (1995) 79. R. Zwanzig, Proc. Natl. Acad. Sci. USA 94 (1997) 148. M.M. Klosek, B.J. Matkowsky, Z. Schuss, Ber. Bunsenges. Phys. Chem. 95 (1991) 331. R. Du, V.S. Pande, A.Y. Grosberg, T. Tanaka, E.I. Shakhnovich, J. Chem. Phys. 108 (1998) 334. J. Shimada, E. Shakhnovich, Proc. Natl. Acad. Sci. USA 99 (2002) 11175. T. Veitshans, D.K. Klimov, D. Thirumalai, Folding Des. 2 (1997) 1. D.K. Klimov, M.R. Betancourt, D. Thirumalai, Folding Des. 3 (1998) 481. E. Shakhnovich, A.M. Gutin, Protein Eng. 6 (1993) 793. E. Shakhnovich, Phys. Rev. Lett. 72 (1994) 3907. D.K. Klimov, D. Thirumalai, Proc. Natl. Acad. Sci. USA 97 (2000) 2544. J.N. Onuchic, N.D. Socci, Z. Luthey-Schulten, P.G. Wolynes, Folding Des. 1 (1996) 441. S.E. Jackson, A.R. Fersht, Biochemistry 30 (1991) 10428. D. Thirumalai, D.K. Klimov, Folding Des. 3 (1999) R112.
258
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
[29] L. Li, L.A. Mirny, E.I. Shakhnovich, Nat. Struct. Biol. 7 (2000) 336. [30] D.E. Otzen, A.R. Fersht, Biochemistry 37 (1998) 8139. [31] C. Clementi, P.A. Jennings, J.N. Onuchic, J. Mol. Biol. 311 (2001) 879. [32] A.R Viguera, L. Serrano, M. Wilmanns, Nat. Struct. Biol. 3 (1996) 874.
[33] L. Li, E. Shakhnivich, J. Mol. Biol. 306 (2001) 121. [34] D. Perl, C. Welker, T. Schindler, K. Schroder, M.A. Marahiel, R. Jaenicke, F.X. Schmid, Nat. Struct. Biol. 5 (1996) 229. [35] R. Koradi, M. Billeter, K. Wuthrich, J. Mol. Graphics 14 (1996) 51.
Progressing from folding trajectories to transition state ensemble in proteins D.K. Klimov a
a,*
, D. Thirumalai
b
Bioinformatics and Computational Biology Program, School of Computational Sciences, George Mason University, Manassas, VA 20110, USA b Department of Chemistry and Biochemistry, Institute for Physical Science and Technology, University of Maryland, Bldg 085, College Park, MD 20742, USA Received 5 March 2004; accepted 25 June 2004 Available online 27 August 2004
Abstract A general procedure for determining the nature of the transition state ensemble (TSE) in proteins is used to probe the heterogeneity of TSE in two models of a/b proteins. The sequences, which differ only in the chain connectivity, are modeled using a-carbon representation. The structures in the TSE are determined using a progress variable clustering (PVC) algorithm that does not depend on the choice of the reaction coordinate. Folding simulations starting from these transition state structures partition with nearly equal probability into the native state or unfolded state. Our fully automated PVC algorithm, which utilizes only the time series in the folding trajectories, is also computationally efficient. The results for the model a/b proteins show that the TSE can be divided into a small number of clusters (two or three) consisting of structurally related conformations. For both sequences the conformations in the dominant transition cluster are compact and native-like with the secondary structure elements fully formed. Although none of tertiary contacts is completely formed in the transition region, few occur with high probability. This suggests that there are multiple (more than one) nuclei that connect the unfolded and native states. Overstabilizing the high probability nucleation contacts retards folding rates. The results for both the models also suggest that, if the transition region is close to the native state, then permutation of the secondary structural elements does not alter the nature of the transition state. Ó 2004 Elsevier B.V. All rights reserved. Keywords: Protein folding; Transition state; Multiple folding nuclei; Progress variable clustering method
1. Introduction The transition from unfolded conformations to the folded state of proteins and RNA involves multiple pathways, in which a number of non-covalent interactions form and break at early times. For two state proteins, near transition state (TS) region a fraction of relatively stable (‘‘nucleating’’) contacts form. Upon formation of the nucleating contacts, which is a mixture of both long and short range contacts (range here refers to the distance along the sequence), the transition to the *
Corresponding author. Tel.: +3014054786; fax: +3013149404. E-mail address: [email protected] (D.K. Klimov).
0301-0104/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.chemphys.2004.06.071
native state occurs with overwhelming probability [1]. The consolidation of the entire structure immediately after the nucleating event is accompanied by the global chain compaction, so that folding and collapse are nearly synchronous. Because of the complexity of protein folding [2–8] it is difficult to obtain a simple reaction coordinate to describe the kinetics of protein folding. Given the multidimensional nature of the free energy landscape physical arguments are often used to obtain a low dimensional pseudo-reaction coordinates [9–11]. For example, for simplified protein models a fraction of native contacts, Q, which measures the assembly of native interactions, has been suggested as a suitable reaction coordinate [9,12]. Several studies demonstrated
252
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
the utility of Q in estimating the free energy barriers for folding as well as the approximate location of the transition state for small two-state folding proteins [9]. In general, it has been argued that Q is an appropriate reaction coordinate provided the protein sequence is highly optimized to fold to the native structure. Unlike in simple chemical reactions involving bond breakage or simple conformational transitions (boat to chair transition in cyclohexane, for example), the very definition of the TS (more precisely, the transition region) is ambiguous. Theoretical arguments suggest that, typically, the transition region is broad with superimposed structure so that it is fruitful to visualize a transition state ensemble (TSE) rather a unique TS for protein folding [13,14]. Such a plasticity in the TSE also rationalizes the ability of many proteins to tolerate large number of mutations without compromising the native fold. Experimental determination of the TS is based on /-value analysis pioneered by Fersht [15,16]. Although the interpretation of /-values is not always straightforward, those close unity indicate that a particular residue is highly structured in the TS. Combined experimental and simulation efforts have given considerable insight into the nature of TSE in a number of proteins. Nevertheless, the straightforward implementation of /-value analysis has been of limited use in the simulations of protein models including the all-atom ones due to considerable computational efforts in reliable determination of thermodynamic and kinetic properties of wild-type proteins [11]. Moreover, the determination of free energies requires the use of suitable order parameter, which can be ambiguous [17]. The computational challenge is to establish the nature of the folding TSE using only the time series generated by suitable dynamics. In other words, given the dynamics of all the particles starting from initially denatured state till the first passage time is reached for a large number of events how can one determine the TSE without making any assumption about the underlying reaction coordinate? Because all-atom molecular dynamics simulations cannot be used to generate a single trajectory for any protein till it reaches the native state starting from unfolded conformation, most of the analysis of TSE has been done using coarse-grained models. In our previous studies we developed a method based on clustering the structures in trajectories and computing TSE by using the ‘‘progress variable’’ d = t/s1i, where s1i is the first passage time for the trajectory i and t s1i [10,14]. Because this method relies both on the progress variable d and clustering of structures, we refer to it as PVC method. In implementing PVC method, we defined the TS structure for a given trajectory as the one which contains the minimal set of stable native contacts, the formation of which ensures rapid passage to the native structure. In practice, we seek a collection of conformations, a putative TSE, which corresponds to sharp rise in
the fraction of native contacts as a function of d. By computing these structures for a large number of independent folding transitions, a meaningful description of the folding TSE may be obtained. An alternative definition of TSE follows from the work of Klosek et al. [18]. In this approach the TS region is taken to be the locus of structures (stochastic separatrix) that have equal probability to lead either to unfolded U or native N states, i.e., the fraction of folding (‘‘forward’’) trajectories pfold 0.5. This idea was first numerically implemented by Du et al. [19] to describe folding kinetics in lattice models without side chains. It has been subsequently applied to more complex lattice models with side chains and other systems [10,20]. Although neither the PVC method nor the pfold technique relies on the choice of a reaction coordinate, they are both subject to certain limitations. The first provides approximate characteristics of the TSE, because its position along the folding pathway is determined from the average folding trajectory. The PVC method is computationally efficient, because it utilizes only the folding trajectories without the need for additional simulations to compute TSE. The stochastic separatrix (pfold) method assumes that at T TF a protein, upon crossing TS, subsequently folds to the native state or unfolds with approximately equal probability [19]. Therefore, if one initiates folding simulations with TSE conformations, the fraction of molecules pfold that finds the native state prior to unfolding should be 0.5. Obvious advantage of pfold method is that the saddle point on the free energy landscape is located by using a physically reasonable (but by no means unique) definition of TS. Yet the obvious computational inefficiency of the pfold method severely limits its practical implementation. In addition, there is some ambiguity in characterizing the unfolded basin of attraction that can compromise the numerical implementation. The purpose of this paper is to compare the TSEs determined by both methods for the coarse-grained models with a/b protein motif. Both coarse-grained Ca-models of a/b proteins differ only in the connectivity of the secondary structural elements (SSEs). For both of those we find that the two methods that use completely different definitions of the TSE give a consistent picture. Although there are appreciable differences in the heterogeneity of folding pathways in the two models the nature of the TSE are quite similar. This surprising result is rationalized in terms of the nearness of the TSE to the native state. 2. Methods 2.1. A coarse-grained model and simulation details for a/b-protein We used a coarse-grained off-lattice Ca-carbon model for an a/b protein. The model sequence consists of
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
N = 38 connected beads of three types, hydrophobic B, hydrophilic L, and neutral N. The potential energy includes contributions from bond-length and bond-angle potentials, the dihedral angle potential, and non-bonded potential described in detail elsewhere [21,22]. The average separation between Ca-carbons along the sequence is ˚ . A sequence is determined by the connectivity a = 3.8 A of residues and the set of random factors k in the nonbonded potential [21]. Two sequences, A1 and A2, were designed to adopt a/b-motifs in their native state, which contains two b-strands (s1 and s2) and a-helix (h) connected by turns (Fig. 1). The a/b topology is stabilized by proper dihedral angle potentials and distribution of hydrophobic residues and k factors achieved using Monte-Carlo sequence design algorithm [23,24]. We used Langevin dynamics (LD) to simulate the folding of the a/b protein. The thermodynamics of folding is computed in the underdamped limit by applying the velocity form of Verlet algorithm [21]. Multiple histogram method was used to compute several probes, such as the fraction of native contacts Q, the fractions of native contacts between SSEs Qij (i,j = s1, s2, h), and the radius of gyration Rg as a function of temperature T. The folding kinetics is obtained from LD simulations in the overdamped limit (approximately at water viscosity) at the temperature Ts . TF, where TF is the folding temperature. The folding transition temperature is determined using standard methods [10]. The mean first passage time sfp is computed by averaging over 400 trajectories. 2.2. Analysis of folding TSE The structures of the TSE are determined using the clustering (pattern recognition) algorithm [10,13,14]. The conformations from folding trajectories generated at Ts are recorded over regular intervals ds = 120 ps(1 + [s1i/600 ns]), where the first passage time s1i expressed in ns corresponds to the first instance of assembly of all native contacts. We tested different values of ds
253
to ensure that qualitatively similar results are obtained. In total, we used 100 folding trajectories. The recorded conformation k is represented by a vector Ck with the components given by the distances between the residues forming native contacts rmn. The vectors Ck (k = 1,. . .,Ns,i, where Ns,i is the number of conformational snapshots for the trajectory i) are clustered to filter out fast structural oscillations [14]. Clustering procedure for a trajectory i requires choosing the values of two parameters, the cluster radius Rc,i and the tolerance for native contacts Dc,i, which defines the maximum deviation of rmn from the native value r0mn . Four criteria are used to select Rc,i and Dc,i: (i) the ‘‘native’’ cluster, which includes the native conformation, should have all native contacts formed; (ii) there shall be no other clusters beside native, in which all native contacts are formed; (iii) the number of structures in the native cluster must be small ([3–4); (iv) the native cluster must first appear close to s1i (e.g., d J 0.99). The optimal values of Rc,i and Dc,i are localized in the upper right corner of the plot in Fig. 2. This choice for Rc,i and Dc,i ensures that the number of clusters Nc,i satisfies the condition 1 Nc,i Ns,i, making the clustering procedure meaningful. The stability of contacts until s1i is determined with a certain tolerance, which permits short-lived disruptions. Following our earlier studies [14,25], we assume that the ‘‘nucleation contacts’’ at a given d constitute the minimal set of native contacts that remain stable until s1i. To establish the location of the putative TSE region, we calculated the kinetic probability Pq(d) averaged over all folding trajectories that a contact q is part of the minimal set of stable native contacts at a given d. From this, we obtain the average (over all native contacts) probability P(d) for a native contact to be part of a minimal set of stable contacts. The value of dTS, at which dP/ dd begins steep rise, is taken to be the TS region. The entire procedure of locating TSE, starting from the clustering of folding trajectories to computing P(d), does not depend on the protein model and is fully automated.
Fig. 1. Ribbon diagrams of native conformations for sequences A1 and A2 visualized using MolMol [35]. Although both native conformations for A1 and A2 have an a/b motif, their topologies are different. A1 adopts strand–strand–helix (b–b–a) fold, whereas the native state of A2 is strand– helix–strand (b–a–b).
254
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
Fig. 2. Diagram of clustering parameters values (Rc,i,Dc,i) meeting the four requirements for cluster procedure described in the text (shaded area). The optimal values of Rc,i and Dc,i correspond to the upper right corner of the shaded area. Efficient elimination of fast structural fluctuation is achieved when Rc,i is as large as possible. The data are plotted for a generic trajectory for sequence A2.
Thus, our method uses only the time series for a number of folding trajectories to determine TSE.
3. Results 3.1. Folding transition state ensemble The TS region has been located using the kinetic probability P(d) (see Section 2). The derivative dP/dd (and P(d), data not shown) remains close to zero for the most part of folding trajectories, implying the absence of stable native interactions (inset to Fig. 3(a)). The range of d, corresponding to the beginning of dramatic rise in dP/dd, is associated with the TS region. After tests of various d we concluded that TS in both sequences, A1 and A2, is located at dTS = 0.9 (inset to Fig. 3(a)). The particular value of dTS is not important as long as dTS is kept within a certain interval (0.85 [ dTS [ 0.95). However, at d < 0.85 a large portion of TS conformation is completely unstructured, whereas at larger d > 0.95 the entire TSE is compact and homogeneous. Propagation of stable native contacts Pq(d) for both sequences (Fig. 4) reveals that, in contrast to rapid assembly of local native helical contacts, stable interactions between SSEs emerge only relatively late in folding process (large d) following the formation of compact transient intermediates. Furthermore, stable tertiary contacts across the sequences form at roughly similar d. As a result the TS is compact and the folding nuclei are dispersed across the sequence [6,26]. These findings imply that the formation of stable native folded structure requires tertiary interactions [27]. Similar propaga-
tion of stable native interactions in both sequences implies that the differences in topological connectivity are not evident in the TSE. This results is a consequence of convergence of folding pathways prior to crossing the TSE. To get a better insight into the formation of native interactions in TSE we consider the probability Pn that the residue n is structured in the TSE (i.e., is engaged in nucleation contacts). From the profiles of Pn at several d, including dTS = 0.90 (Fig. 3) the following conclusions can be drawn: (a) The majority of residues in A1 and A2 are involved to a varying degree in the folding nuclei (with the exception of a few, mostly in turn regions, which do not make native interactions), but none of the contacts is fully structured (Pn = 1) [14,28,29]. As a result the nucleus is delocalized and the TS structures are compact [6,26]. (b) Although most of residues are involved in the folding nuclei to some extent, there are regions, which are more structured than the others. Due to topological connectivity, residues involved in nucleation are often located near the turns, connecting the SSE as predicted earlier [13]. This is seen for the residues in the A1 hairpin adjacent to the s1–s2 turn. In A2, the residues in the strand s1 and helix are more structured than those in s2, which suggests that the TS in A2 is more ‘‘polarized’’ than in A1. The interactions formed in the TS are expected to strongly affect folding times. To verify this expectation we selected in A1 eight TS hydrophobic contacts with high Pq(dTS). Four mutants were designed by increasing
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
(a)
(b) Fig. 3. Kinetic probabilities Pn for the residue n at three values of d close to the TS value dTS = 0.90, for A1 (a) and A2 (b). Pn, which estimates the degree to which a particular residue n is structured in the TS, shows that folding nuclei are diffused as the majority of residues participate to a varying degree in the TSE. The inset in (a) shows dP/ dd for A1 (see Section 2). The value of dTS = 0.90 marked by dashed vertical line corresponds to the beginning of dP/dd growth, which is associated with crossing TS. dP/dd for A2 shows similar behavior. Data are calculated using 100 folding trajectories.
the strength of these interactions by 0.63, 1.25, 1.88 and 2.50 kcal/mol (variation by about 8–30%). Compared to the WT the folding times of the first two mutants decrease (to 0.93sF and 0.89sF). However, further stabilization of TS interactions (last two mutants) impedes folding and the folding time increases to 1.05sF and 1.27sF. Therefore, the strength of TS interactions must be optimal as too strong TS contacts may lead to entanglement traps, which in turn slow down folding [14]. 3.2. Structural diversity of the TSE To assess structural diversity of TSE, we applied the clustering analysis to TSE structures [25]. The A1 TSE
255
can be divided into three distinct clusters. The first, TS1, is compact and native-like. The ratio of the radius of gyration of TS1 conformations to its native value Rg =RN g is 1.0. TS1 contains, on an average, about twothirds of native contacts (Q = 0.68). The secondary structures in TS1 are completely formed. Most of the tertiary interactions between SSE are formed with the hairpin s1–s2 being the most structured region. TS1 appears in 70% of all folding trajectories. The structures in the TS2 cluster are less compact ðRg =RN g 1:3Þ and structured than TS1. In TS2 the hairpin is mostly formed, while there are no significant interactions between the helix and strands. The secondary structure in TS2 is not as fully formed as in TS1 (the b-strand and a-helix contents are between 0.5 and 0.6). About 25% of folding trajectories pass through TS2. The structure of rare TS3 is similar to TS2, but differs in the location of the structured region. In TS3, the tertiary interactions are formed only between the strand s2 and h. There are two dominant TS clusters in A2, the first of which, TS1, is very similar to TS1 of A1. This TS is compact ðRg =RN g ¼ 1:0Þ and native-like (Q = 0.60). About 80% of unfolded conformations reach the native state through TS1, while the remaining 20% of molecules fold via the less structured TS2 ðQ ¼ 0:22 and Rg =RN g ¼ 1:4Þ. More importantly, the tertiary interactions between SSEs in TS2 are found only between the strand s1 and helix. 3.3. Application of pfold to TSE For two-state folders the TSE may be defined as a set of structures, from which the probabilities to reach native N or unfolded U states are approximately equal. Folding trajectories initiated with TSE conformations should partition into those, which reach N prior to unfolding (‘‘forward’’ trajectories) and those, which initially unfold. The fraction of ‘‘forward’’ trajectories, pfold, should be approximately 0.5. This requirement is physically reasonable operational definition of the TSE [19]. To assess whether the TSE identified by the PVC method is consistent with the pfold definition, we calculated pfold for the TS structures obtained at dTS = 0.90. In all, we recorded 129 putative TS conformations for A1 and 152 for A2 and selected every fifth conformation. Each such conformation serves as a starting point for M = 100 folding trajectories generated at T = Ts. Each trajectory either reaches the native state (Q = 1.0) or leads to U (i.e., when the number of tertiary contacts Cu < 5). The value of Cu was varied to ensure that the resulting pfold values do not change significantly. For the kth putative TS conformation we calculated the coefficient pkfold ¼ M kF =M, where M kF is the number of ‘‘forward’’ trajectories. Then pfold is taken to be an
256
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
(a)
(b) Fig. 4. The kinetic probabilities Pq as a function of d averaged over 100 trajectories for A1 (a) and A2 (b). Pq gives the probability that a native contact q belongs to a minimal set of stable native contacts at a given d which survives until the first passage time. The color codes to Pq are given on the right. The plots show a remarkable heterogeneity in the formation of native contacts. Helical q become stable very early in the folding process, while all other native contacts irrespective of their location stabilize close to reaching the native state, i.e., the TS occurs late in the folding process. Surprisingly, plots for A1 and A2 are qualitatively similar. The native contact indexes q are arranged according to the type of secondary structure and sequence separation.
average over all pkfold . The values of pfold for A1 and A2 are equal to 0.51 and 0.59, respectively. Thus, the descriptions of the TSE obtained by the stochastic separatrix method and the PVC method are consistent.
3.4. Chain connectivity and the nature of TSE The dominant transition state (TS1) for both the models have similar characteristics (Fig. 1) despite the fact
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
that A1 and A2 have different connectivity between the SSEs. Experimentally the variations in the chain connectivities can be achieved by circular permutation. Whether or not the TSE shifts upon a circular permutation seems to depend on the protein topology and the precise location of permutation. It has been suggested that the nature of the TSE does not change upon altering chain connectivity for CI2 [30,31], but does so dramatically for the SH3-domain of a-spectrin [32,33]. Our simulations suggest that, if the structures in the dominant TS resemble closely that of the native state, then circular permutation would have little effect on the TSE. In both the model sequences structures in the dominant TSs are nearly nativelike with fully formed secondary structures. Such a highly native-like structured TS is only possible, if earlier events lead to specific collapse. The similarity in the structural details of the dominant TSs of two models suggests that the TSE is crossed in the last stages of folding. Experimentally this would imply that if the Tanford b values are close to unity, as is the case in a class of cold-shock proteins [34], then the TSE would not alter appreciably upon permutation provided that circular permutation does not alter the fold itself. For this class of proteins we predict that the TSE would be relatively homogeneous for all the mutants.
4. Conclusions In this study we have used two independent methods for analyzing folding TSE. The stochastic separatrix (pfold) method, which provides a physically motivated description of TSE, confirms the location of TS established using PVC method. Through the analysis of hundreds of folding trajectories we show that dominant folding TS exists for both sequences, A1 and A2. Most molecules (70–80%) fold through this TS, which is compact and topologically resembles the native fold. In the dominant TS the SSEs are fully formed. This conclusion is consistent with the calculated time scales for collapse sc and folding sF. The ratio sF/sc is 2.7 for A1 and 2.4 for A2. These results together with the TS analysis show that the dominant folding TS (corresponding to the ratelimiting step in folding) is native-like, compact, and involves interactions (to a varying degree) between all SSEs. This picture of TSE resembles the notion of delocalized folding nuclei diffused across the native structure [6,9,26]. Because the rate-limiting step occurs late in the compact phase, the dominant folding TSs for A1 and A2 are quite similar, even though their native states are topologically different. Despite the existence of a dominant single TS, there is still certain diversity in the TSE as 20– 30% of molecules fold through alternative TSs. The structure of TS2 is polarized, because native interactions are predominantly formed near the N-terminal (the bhairpin in A1 or s1–h construct in A2).
257
This study shows that a combination of two methods based on d progress variable (a PVC method) and the notion of stochastic separatrix (pfold method) offers a powerful tool for locating and analyzing folding TSE. Because both the methods are general and do not rely on a choice of a folding reaction coordinate, they may be applied to variety of protein models, from simple lattice to all-atom representations. The computationally more efficient PVC method is easier to use.
Acknowledgment This work was supported in part by the National Science Foundation Grant No. NSF CHE-0209340.
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]
[16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]
P.G. Wolynes, Proc. Natl. Acad. Sci. USA 93 (1996) 14249. A.R. Fersht, Proc. Natl. Acad. Sci. USA 92 (1995) 10869. Z. Guo, D. Thirumalai, Biopolymers 36 (1995) 83. E. Shakhnovich, V. Abkevich, O. Ptitsyn, Nature 379 (1996) 96. A.R. Fersht, Curr. Opin. Struct. Biol. 7 (1997) 3. P.G. Wolynes, Proc. Natl. Acad. Sci. USA 94 (1997) 6170. J.N. Onuchic, Z.A. Luthey-Schulten, P.G. Wolynes, Ann. Rev. Phys. Chem. 48 (1997) 545. L. Mirny, E. Shakhnovich, Ann. Rev. Biophys. Biomol. Struct. 30 (2001) 361. N.D. Socci, J.N. Onuchic, P.G. Wolynes, J. Chem. Phys. 104 (1996) 5860. D.K. Klimov, D. Thirumalai, Prot. Struct. Funct. Genet. 43 (2001) 465. C. Clementi, H. Nymeyer, J.N. Onuchic, J. Mol. Biol. 298 (2000) 937. C. Clementi, P. Jennings, J. Onuchic, Proc. Natl. Acad. Sci. USA 97 (2000) 5871. Z. Guo, D. Thirumalai, Folding Des. 2 (1997) 377. D.K. Klimov, D. Thirumalai, J. Mol. Biol. 282 (1998) 471. A.R. Fersht, Structure and Mechanism in Protein Science: A Guide to Enzyme Catalysis and Protein Folding, Freeman, New York, 1998. A.R. Fersht, Curr. Opin. Struct. Biol. 5 (1995) 79. R. Zwanzig, Proc. Natl. Acad. Sci. USA 94 (1997) 148. M.M. Klosek, B.J. Matkowsky, Z. Schuss, Ber. Bunsenges. Phys. Chem. 95 (1991) 331. R. Du, V.S. Pande, A.Y. Grosberg, T. Tanaka, E.I. Shakhnovich, J. Chem. Phys. 108 (1998) 334. J. Shimada, E. Shakhnovich, Proc. Natl. Acad. Sci. USA 99 (2002) 11175. T. Veitshans, D.K. Klimov, D. Thirumalai, Folding Des. 2 (1997) 1. D.K. Klimov, M.R. Betancourt, D. Thirumalai, Folding Des. 3 (1998) 481. E. Shakhnovich, A.M. Gutin, Protein Eng. 6 (1993) 793. E. Shakhnovich, Phys. Rev. Lett. 72 (1994) 3907. D.K. Klimov, D. Thirumalai, Proc. Natl. Acad. Sci. USA 97 (2000) 2544. J.N. Onuchic, N.D. Socci, Z. Luthey-Schulten, P.G. Wolynes, Folding Des. 1 (1996) 441. S.E. Jackson, A.R. Fersht, Biochemistry 30 (1991) 10428. D. Thirumalai, D.K. Klimov, Folding Des. 3 (1999) R112.
258
D.K. Klimov, D. Thirumalai / Chemical Physics 307 (2004) 251–258
[29] L. Li, L.A. Mirny, E.I. Shakhnovich, Nat. Struct. Biol. 7 (2000) 336. [30] D.E. Otzen, A.R. Fersht, Biochemistry 37 (1998) 8139. [31] C. Clementi, P.A. Jennings, J.N. Onuchic, J. Mol. Biol. 311 (2001) 879. [32] A.R Viguera, L. Serrano, M. Wilmanns, Nat. Struct. Biol. 3 (1996) 874.
[33] L. Li, E. Shakhnivich, J. Mol. Biol. 306 (2001) 121. [34] D. Perl, C. Welker, T. Schindler, K. Schroder, M.A. Marahiel, R. Jaenicke, F.X. Schmid, Nat. Struct. Biol. 5 (1996) 229. [35] R. Koradi, M. Billeter, K. Wuthrich, J. Mol. Graphics 14 (1996) 51.