INTRODUCTION

TOP ABSTRACT INTRODUCTION MODELS AND SIMULATION METHOD RESULTS DISCUSSION REFERENCES

In our laboratory we are using computer simulation to understand the basic physical principles of protein aggregation, a cause or associated symptom of a number of lethal disorders including Alzheimer's, Parkinson's, and prion diseases (Clark and Steele, 1992; Eaton and Hofrichter, 1990; Gallo et al., 1996; Massry and Glasscock, 1983; Moore and Melton, 1997; Selkoe, 1991; Simmons et al., 1994). Our long-term goal is to learn if there are certain universal molecular-level mechanisms responsible for protein aggregation and fibril formation. Because the association of -strands, or the conversion of -helices to -strands often precedes or accompanies protein aggregation (Benzinger et al., 1998, 2000; Burkoth et al., 1998; Esler et al., 1996, 2000; Lazo and Cowning, 1998; Lynn and Meredith, 2000; Sunde et al., 1997; Zhang et al., 2000), we have focused part of our efforts on the development of models of -strand proteins. These models must be capable of capturing the essential physical features of real -strand proteins and at the same time be simple enough to allow the simulation of multi-protein systems with current computer capability.

Low-resolution or simplified protein models are best suited for the construction of a -strand model to be used in simulations of protein folding because they allow us to study the folding behavior over relatively long time scales. These models provide numerous insights into the thermodynamic and kinetic properties of protein folding (Chan and Dill, 1994; Dokholyan et al., 1998, 2000; Guo and Brooks, 1997; Guo and Thirumalai, 1995, 1996; Gupta and Hall, 1997, 1998; Kolinski et al., 1995, 1999; Lau and Dill, 1989; Miller et al., 1992; Nymeyer et al., 1998; Pande and Rokhsar, 1998; Shea et al., 2000; Skolnick and Kolinski, 1991; Zhou and Karplus, 1997a,b, 1999). We recently introduced three minimalist models of four-strand antiparallel -strand peptides: the -sheet, the -clip, and the -twist (Jang et al., 2002). Discontinuous molecular dynamics (DMD) simulations (Alder and Wainwright, 1959; Rapaport, 1978; Smith et al., 1996) on these three models were performed to determine how the thermodynamic properties of an isolated peptide vary with temperature. Despite these models' simplicity, they undergo a complex set of protein transitions similar to those observed in experimental studies on real proteins (Ptitsyn, 1995). Starting from high temperature, these transitions include a collapse transition, a disordered-to-ordered globule transition, a folding transition, and a liquid-to-solid transition. The thermodynamics results presented in our previous paper and the kinetics results presented here set the stage for a later paper in which we will examine the folding and assembly of several of these model -strand proteins into fibrils.

The folding of -strand proteins is of current interest not only because isolated -strand peptides tend to adopt a -sheet fibril conformation when they aggregate, but also because -strands are often secondary structure elements in proteins. Unfortunately, unlike the formation of -helices, the formation of even isolated -strand peptides is not easy to study experimentally due to difficulties associated with the propensities of the largely insoluble -structures to precipitate. This has prompted the development of designed -strand peptides that are more amenable to fundamental protein folding investigations. Experimental studies of such peptides include analyses of the structure and stability for de novo designed two-strand -sheets known as -hairpin (Blanco et al., 1998; Ramírez-Alvarado et al., 1996), three-strand -sheets known as betanova (de Alba et al., 1999; Koepf et al., 1999; Kortemme et al., 1998), and four-strand -sheets known as betabellin (Lim et al., 2000). Computer simulation studies have been conducted on the thermodynamics and kinetics of formation of 16-residue -hairpins (Dinner et al., 1999; García and Sanbonmatsu, 2001; Lee and Shin, 2001; Pande and Rokhsar, 1999; Kolinski et al., 1999; Zagrovic et al., 2001; Zhou and Linhananta, 2002; Zhou et al., 2001) taken from the C-terminal fragment (41-56) of the B1 domain of protein G (Muñoz et al., 1997), the folding kinetics of betanova (Bursulaya and Brooks, 1999; Ferrara and Caflisch, 2000), the folding thermodynamics and kinetics of a 46-mer -barrel (Guo and Brooks, 1997; Guo and Thirumalai, 1995), and the dynamics of a six-stranded lattice model of Greek key -barrel (Kolinski et al., 1995; Skolnick et al., 1989).

In this paper we investigate the folding kinetics of the three different four-strand antiparallel -strand peptides: the -sheet, the -clip, and the -twist, whose thermodynamic properties were studied previously by us (Jang et al., 2002). These -strand peptides each have 39 connected residues (beads) and different native state conformations as shown in Fig. 1. Nonbonded beads can interact through a hybrid Go-type potential (Go and Taketomi, 1978, 1979; Taketomi et al., 1975; Ueda et al., 1978) modeled as a square-well or square-shoulder potential depending on the value of the bias gap parameter g. DMD simulations were performed on these three model systems at intermediate values of the bias gaps ranging from 0.7 to 1.1. The progress parameters, the squared radius of gyration R, the fraction of native contacts within the protein as a whole Q, and between specific strands Q_ab were monitored over the course of the kinetic simulations. The bias gap measures the difference in interaction strength between the native and non-native contacts; the larger it is, the more the native state is favored over the non-native state. Intermediate values of the bias gap are thought to be the most representative of real proteins in their equilibrium and dynamic behavior. The model protein's bias gap is an artificial measure of its preference for the native state. By exploring how variations in the bias gap influence the types of dynamic phase behavior observed, we can get a feeling for how a real protein's preference for its native state, as measured for example by the energy difference between the denatured and native state, is manifested in the protein's dynamic behavior and vice versa.

View larger version (36K): [in this window] [in a new window]   FIGURE 1   Stick (top) and topology (bottom) diagrams of the global energy minimum structures: (A) the -sheet; (B) the -clip; and (C) the -twist model proteins.

Highlights of our results are the following. We find that the -sheet exhibits a fast-track folding pathway without becoming trapped in any intermediates for all the bias gaps investigated. Its kinetic folding mechanism can be described by the diffusion collision model (Karplus and Weaver, 1979, 1994). In contrast, the -clip and -twist exhibit multiple folding pathways that include trapping in intermediates and direct folding to the native state depending on the size of the bias gap. For the -clip and -twist models at g = 0.7, two intermediates I₁ and I₂ are observed in the folding pathways. The intermediate I₁ is a partially ordered globule, and the intermediate I₂ is an ordered globule that can be regarded as a molten globule as described in the thermodynamic study. The intermediate I₁ can fold either directly into the native state or via the intermediate I₂, which then folds into the native state by elongating its strands. No fast-track folding pathway is observed. The kinetic folding mechanism for the -clip and -twist at g = 0.7 can be described by the hydrophobic-collapse theory (Dill et al., 1995). For the -clip and -twist models at g = 0.9, the kinetic folding trajectories include the fast-track folding pathway and trapping in an intermediate I. The intermediate I is an ordered globule as was seen in the g = 0.7 model result. The random coil can fold directly to the native state or via the intermediate I. We expect that as g increases, more of the trajectories will exhibit the fast-track kinetic folding pathway. The kinetic folding mechanism for the -clip and -twist at g = 0.9 can be described by a combination of the diffusion collision model and the hydrophobic-collapse theory. The folding speed of the model proteins with different native state topologies strongly depends on the contact order in the native state (Plaxco et al., 1998). Models with a high topological complexity exhibit slower folding to the native state than models with a topologically simple structure. Moreover, the former is prone to kinetic intermediates.

In the following section, a full description of the three models and the simulation method are presented. Next is the results section, which is divided into three subsections; the first presents the results for the folding kinetics, the second presents a description of intermediates and folding pathways, and the third discusses how the folding speed depends on the topology of the model proteins. The paper concludes with a discussion of the key findings in the last section.

	MODELS AND SIMULATION METHOD

TOP ABSTRACT INTRODUCTION MODELS AND SIMULATION METHOD RESULTS DISCUSSION REFERENCES

We consider three different off-lattice protein models whose native states are four-strand antiparallel -strand peptides. The global energy minimum structures for the three different -strand models are shown in Fig. 1. Stick diagrams for the -sheet, -clip, and -twist are shown at the top of the figure. Topology diagrams representing a top view of the three proteins are shown at the bottom of the figure. In the topology diagrams, thick arrows indicate chain connectivity at the top of the strands, and thin arrows indicate chain connectivity at the bottom of the strands. It is immediately apparent that in the native state, the -sheet is two-dimensional (planar), whereas the -clip and -twist have well-defined three-dimensional four-barrel structures. The total number of native contacts, N, for each model chain and the total number of inter-strand native contacts, N, for each model where a and b are the strand number, are summarized in Table 1. Note that N for the -sheet equals the sum of all N, but N for the -clip and -twist does not, because native contacts on turn residues are not inter-strand native contacts.

Each model protein contains 39 connected beads, each representing an amino acid residue that can be regarded as being localized at the C atom. The details of the model potential and the bias gap g are described in an earlier publication (Jang et al., 2002).

The folding kinetics of the three different -strand peptides are studied here using the DMD algorithm (Alder and Wainwright, 1959; Rapaport, 1978; Smith et al., 1996). Proteins are quenched from the denatured or random coil state equilibrated at high temperature to the temperature of interest. After quenching, the proteins undergo a conformational change toward the energy minimum state, i.e., the native state. We investigated the time series of folding behavior from the random coil state at T* = 2.0 to the native state at T* = 0.2, 0.24, and 0.28 for the -sheet, the -clip, and the -twist, respectively. The temperature used in this paper is a dimensionless reduced temperature, T* = k_BT/, where k_B is Bolzmann's constant, T is the absolute temperature, and  is the energy parameter. Our thermodynamic study on the same models (Jang et al., 2002) shows that at T* = 2.0, the three -strand peptides are completely random coils. The three target temperatures for each model are selected to ensure that the native state is highly populated at that temperature. The initial configurations for the kinetic simulations were generated from equilibrium simulations at T* = 2.0. The coordinates and velocities of the beads in the random coil state were obtained every 10⁵ collisions during the equilibrium simulation to create independent initial conformations for the quenching studies. This large collision interval ensures that the sampled configurations are independent. Simulations were performed for three different values of the bias gaps, g = 0.7, 0.9, and 1.1, except for the -sheet where the simulations were performed at g = 0.9 and 1.1. For each bias gap model, statistical averages were obtained over 100 independent simulations.

The total time duration for the simulation strongly depends on the topology of the model and the bias gap. The time used in the simulations is expressed in terms of reduced time, t* = t (Zhou and Karplus, 1999), where t is the real time,  is the energy parameter, M is the mass of the bead, and  is the bead diameter. For all three -strands, simulations were performed for ~10⁸ collisions, which is equivalent to t* = 4 × 10⁵ for the -sheet model and t* = 2 × 10⁵ for the -clip and -twist models. This indicates that the simulation speed for the -sheet is twice as fast as that of the other two models. The simulation speed at large bias gaps is slightly faster than at small bias gaps for any given model, indicating that the larger the bias gap, the faster the protein folds. To highlight the folding events that occur early in the simulation, a logarithmic timescale was used to examine early folding behavior of proteins as well as folded behavior in the native state. Simulation results were recorded at linearly spaced time intervals of t* = 0.1 for t* < 10 and at logarithmically spaced time intervals of log t* = 0.01 for t* > 10.

To determine the compactness of the chains during the simulations, the time-dependent squared radius of gyration, R (t) was determined where

(1)

with N equal to the number of beads. The time-dependent bead positions, r_i(t) (i = 1 to N), are obtained by shifting the center of mass coordinates to the origin. The radius of gyration is widely used as an order parameter in studies of the polymer order-disorder collapse transition. At low temperatures, collapsed polymers tend to form a high-density conformation yielding a small value of the radius of gyration. In contrast, the squared radius of gyration for native -strand peptides is larger than that in the collapsed state, because the native -strands are long and have an elongated molecular shape. This indicates that a minimum value of the squared radius of gyration does not imply that the -strands are in the native state. Thus, to characterize the nativeness of -strands as collapsing progress, an alternative measure of the time-dependent squared radius of gyration is considered. This is the normalized squared radius of gyration, F^*_R(t), which is defined by

(2)

where R(eq) is the native value of the squared radius of gyration. The normalized squared radius of gyration measures the departure of the chain's squared radius of gyration from its value in the native state. It is slightly less than 1 when the chain is loosely packed or a random coil (R(t) > R(eq)), equals 0 when the chain is in the native state (R(t) = R(eq)), and is less than 0 when the chain collapses to very compact structure (R(t) < R(eq)). Each -strand peptide has different values of R(eq) depending on the temperature and the size of the bias gap. The native values of R for the three -strands at different bias gaps are obtained from the equilibrium simulations at T* = 0.2, 0.24, and 0.28 for the -sheet, the -clip, and the -twist, respectively. These values are summarized in Table 2 along with those for the global energy minimum structures. The R(eq) values for the native -strands are smaller than the R values in the global energy minimum structures, R(GM), because the -strands are slightly bent or kinked in the native state, especially for the smaller bias gap models. This is due to the attraction between non-native contacts for g < 1, which disturbs the propensity of the -strands to be straight. Note that the ratio of the native R value to that in the global energy minimum structure, R(eq)/R(GM), is much smaller for the -sheet than for the other two models. This indicates that the planar shape of the -sheet can be easily bent due to the non-native attractions, giving the native state a slightly rolled molecular structure. In fact, most real -sheet structures are twisted as opposed to being perfectly planar or elongated molecular structures (Clothia, 1973; Pauling and Corey, 1951).

The progression of a conformation toward the native state can be monitored by introducing the time-dependent fraction of native contacts formed (Lazaridis and Karplus, 1997; ali et al., 1994b), Q(t), defined by

(3)

and further breaking down Q(t) into the time-dependent fraction of inter-strand native contacts formed between strands a and b, Q_ab(t), defined by

(4)

where N_native represents the number of native contacts for the whole chain and N is the number of inter-strand native contacts between strands a and b in a given conformation. The fraction of native contacts Q has a value between 0 and 1. When Q = 1 the model system can be regarded as being in the native structure, whereas when Q  0 the chain becomes a random coil that signifies the denatured state.

	RESULTS

TOP ABSTRACT INTRODUCTION MODELS AND SIMULATION METHOD RESULTS DISCUSSION REFERENCES

Characteristics of folding kinetics

We performed DMD simulations to investigate the folding kinetics of the three different -strand peptides. The systems were quenched from the denatured state at T* = 2.0 to the native states at T* = 0.2, 0.24, and 0.28 for the -sheet, the -clip, and the -twist, respectively. The time-dependent normalized squared radius of gyration, F^*_R, is shown in Fig. 2 as a function of the reduced time on a logarithmic scale, log(1 + t*): for the -sheet at g = 0.9 and 1.1 (Fig. 2 A); for the -clip at g = 0.7, 0.9, and 1.1 (Fig. 2 B); and for the -twist at g = 0.7, 0.9, and 1.1 (Fig. 2 C). The quantity F^*_R is an average over 100 independent trajectories that start from different initial configurations. Simulation results are presented up to the reduced time of t* = 10⁵, because for t* > 10⁵ the -peptides are very stable and remain in the folded state. The degree of compactness of the system can be monitored by examining how the F^*_R values change as folding progress. Recall that when 0 < F^*_R < 1, the chain has a larger squared radius of gyration, R, than the native value of R(eq) (see Table 2), indicating that the system is loosely packed or a random coil. For F^*_R = 0, the system has the same R value as in the native state, so that the system can be regarded as being in the native state. For F^*_R < 0, the system exhibits a highly collapsed state, yielding a kinked or rolled molecular structure with smaller value of R than in the native state.

View larger version (17K): [in this window] [in a new window]   FIGURE 2   The average values of the normalized squared radius of gyration, F*R(t), as a function of the reduced time in logarithmic scale, log(1 + t*): (A) the -sheet at g = 0.9 and 1.1; (B) the -clip at g = 0.7, 0.9, and 1.1; and (C) the -twist at g = 0.7, 0.9, and 1.1.

The F^*_R(t) versus log(1 + t*) graph in Fig. 2 A for the -sheet shows that F^*_R > 0 initially, that it oscillates around zero when the native state is reached, and that the folding events for the g = 0.9 and g = 1.1 models are qualitatively the same. The initial conformations of the chains begin to collapse at t* = 10 and fold into the native state by t* = 100. The F^*_R(t) versus log(1 + t*) graphs in Fig. 2, B and C, for the -clip and the -twist, respectively, are strikingly different. In contrast to the -sheet, the F^*_R values for the -clip and the -twist rapidly decrease to a negative minimum value as the chain collapses and then gradually increase to zero as the native state is reached. The negative value for F^*_R means that the R value of the collapsed chain is smaller than that of the native state at the given time. This suggests that the folding process for the -clip and -twist goes through two distinct stages. In the first stage, the initial random coil collapses to a compact structure. In the second stage, strand extension takes place, recovering the elongated molecular shape of the -strand peptide. For the -clip in Fig. 2 B, the minimum value of F^*_R decreases as g decreases, indicating that the smaller the value of g, the more compact the collapsed structure is. This is due to the relatively high overall strength of the non-native attractions among the beads in the smaller gap models. The strand extension is much quicker in the larger gap model. Similar behavior is seen in Fig. 2 C for the -twist.

The time-dependent fraction of native contacts, Q, is shown in Fig. 3 as a function of log(1 + t*): for the -sheet at g = 0.9 and 1.1 (Fig. 3 A); for the -clip at g = 0.7, 0.9, and 1.1 (B); and for the -twist at g = 0.7, 0.9, and 1.1 (C). The quantityQ is an average over 100 independent simulations with different independent initial configurations. For the -sheet, no difference in the formation of native contacts between the g = 0.9 and g = 1.1 models is observed. The formation of native contacts and the chain collapse (compare to Fig. 2 A) are concurrent. This indicates that the collapse time and the folding time are identical for the -sheet. The native states for both the bias gap models are located at Q  1. In contrast, for the -clip in Fig. 3 B and for the -twist in Fig. 3 C, 60-70% of the total native contacts (depending on the size of the bias gap) are built up when F^*_R reaches its negative minimum at around t* = 100 as seen in Fig. 2, B and 2 C. The Q values continuously increase thereafter and begin to saturate slowly to the native value during the strand extension process as F^*_R goes to zero. This reflects the result that for the -clip and -twist the collapse time and the folding time are different. For the -clip, the native states are located at Q  0.82, 0.92, and 0.95 for g = 0.7, 0.9, and 1.1, respectively. As g increases, the Q value in the native state increases, and the speed of native contact formation is much faster. For the -twist, the native states are located at Q  0.88, 0.92, and 0.93 for g = 0.7, 0.9, and 1.1, respectively. The reason that the native Q values are all different for the three -strand models with the same bias gap is that the kinetic simulations are performed at different target temperatures.

View larger version (14K): [in this window] [in a new window]   FIGURE 3   The average values of the fraction of native contacts, Q(t), as a function of the reduced time in logarithmic scale, log(1 + t*), for: (A) the -sheet at g = 0.9 and 1.1; (B) the -clip at g = 0.7, 0.9, and 1.1; and (C) the -twist at g = 0.7, 0.9, and 1.1.

More detailed information on the formation of native contacts can be obtained by considering the fraction of inter-strand native contacts, Q_ab, between strands a and b where the strand number runs from 1 to 4. The -sheet model has three possible strand-to-strand native interactions: 1-2, 2-3, and 3-4 (see the topology diagrams in Fig. 1 A). The time-dependent fraction of inter-strand native contacts, Q_ab, for the -sheet model is shown in Fig. 4 as a function of log(1 + t*) for the bias gaps: g = 0.9 (Fig. 4 A) and g = 1.1 (Fig. 4 B). Each value of Q_ab is averaged over 100 independent simulations. For g = 0.9, Q₁₂ and Q₃₄ (the fraction of inter-strand native contacts experienced by the outer strands) increase before Q₂₃ (the fraction of inter-strand native contacts experienced by the inner strands) increases as the system collapses. This indicates that assembly of the native state for the -sheet at g = 0.9 starts with the formation of native contacts between strands 1 and 2 and between strands 3 and 4, followed by the formation of native contacts between strands 2 and 3. This is due to the fact that strands 1 and 4 are the chain's terminals, which are more flexible and easier to rearrange than the other strands. We find that Q₁₂  Q₂₃  1 at t* = 160, and Q₃₄  1 at t* = 200. However, for g = 1.1 all Q_ab values increase concurrently and Q₁₂  Q₂₃  Q₃₄  1 at t* = 100. The native secondary structure formation for the -sheet at g = 1.1 starts at the same time as the chain collapses.

View larger version (14K): [in this window] [in a new window]   FIGURE 4   The average values of the fraction of inter-strand native contacts, Qab(t), where a and b are the strand number, for the -sheet as a function of the reduced time in logarithmic scale, log(1 + t*), for the bias gaps: (A) g = 0.9; (B) g = 1.1.

In marked contrast to the -sheet, the formation of inter-strand native contacts for the -clip is more complex as shown in Fig. 5 A at g = 0.7 and Fig. 5 B at g = 0.9. As shown in Fig. 1 B (top view topology diagram), there are six possible strand-to-strand native interactions for the -clip: Q₁₂, Q₁₄, Q₂₃, and Q₃₄ correspond to the four peripheral strand-to-strand interactions, and Q₁₃ and Q₂₄ correspond to the two diagonal strand-to-strand interactions. For g = 0.7, the formation of -clip native contacts begins to develop at the ends of the chain; i.e., Q₁₂ and Q₃₄ increase first, as was the case in Fig. 4 A for the -sheet at g = 0.9, because the chain's terminals (strands 1 and 4) are very flexible and easy to rearrange. All values of Q_ab increase initially as the chain collapses and saturate at their respective native values of Q_ab, except for Q₁₄ which continuously increases and reaches a value of Q₁₄  0.85 at t* = 10⁵. The native state is located at t*  10⁵, the time at which Q₁₄ starts to oscillate about 0.85. In the native state, the four peripheral values, Q₁₂, Q₁₄, Q₂₃, and Q₃₄ oscillate about their average value of Q  0.82, and the two diagonal values, Q₁₃ and Q₂₄ oscillate about 0.6, which is much smaller. Because the number of native contacts in the four peripheral strand-to-strand interactions is much larger than that in the two diagonal strand-to-strand interactions (see Table 1), this indicates that the larger the number of strand-to-strand native contacts, the faster the formation of native contacts in the strand-to-strand interactions. For g = 0.9, Fig. 5 B shows similar information on the formation of inter-strand native contacts. The speed of the formation of native contacts is faster than in the g = 0.7 case. The native state is located at t*  10⁴ when Q₁₄ reaches its native value of 0.95. Note that in the native state the diagonal strand interactions Q₁₃ and Q₂₄ oscillate at around 0.8, which is significantly closer to the native value for the whole chain of Q  0.92 than occurs in the g = 0.7 case. This is based on the fact that the larger the bias gap, the more the native state is favored over the non-native state. In fact, the large values of Q₁₃, Q₁₄, and Q₂₄ are essential to maintain the -clip structure, because without these strand-to-strand native interactions, the model -clip reduces to the -sheet.

View larger version (20K): [in this window] [in a new window]   FIGURE 5   The average values of the fraction of inter-strand native contacts, Qab(t), where a and b are the strand number, for the -clip as a function of the reduced time in logarithmic scale, log(1 + t*), for two selected values of the bias gap: (A) g = 0.7 and (B) g = 0.9.

The formation of inter-strand native contacts for the -twist is shown in Fig. 6 A at g = 0.7 and Fig. 6 B at g = 0.9. As shown in Fig. 1 C (top view topology diagram), there are five possible strand-to-strand native interactions for the -twist: 1-2, 1-3, 2-3, 2-4, and 3-4. For g = 0.7 in Fig. 6 A, native contacts between strands 1 and 2 and between strands 3 and 4 form rapidly as folding progresses, but native contacts between strands 1 and 3 and between strands 2 and 4 form slowly. Native contacts between strands 2 and 3 form at an intermediate speed. After the initial rapid increase, Q₁₂ and Q₃₄ oscillate about 0.88 in the native state, whereas Q₁₃ and Q₂₄ continuously increase to 0.88. The native state is located at t*  3.2 × 10⁴ when Q₁₃ and Q₂₄ begin to oscillate about 0.88. For g = 0.9 in Fig. 6 B, it is immediately apparent that there are three fast rates in the formation of inter-strand native contacts (Q₁₂, Q₂₃, and Q₃₄) and two slow rates in the formation of inter-strand native contacts (Q₁₃ and Q₂₄). The quantities Q₁₃ and Q₂₄ increase faster than they do in the g = 0.7 model. The native state is located at t*  10⁴ when Q₁₃ and Q₂₄ oscillate about 0.92. For the -twist, the large values of Q₁₃ and Q₂₄ are key factors in the conformation of the native -twist structure, because without these strand-to-strand native interactions, the model -twist reduces to the -sheet.

View larger version (19K): [in this window] [in a new window]   FIGURE 6   The average values of the fraction of inter-strand native contacts, Qab(t), where a and b are the strand number, for the -twist as a function of the reduced time in logarithmic scale, log(1 + t*), for two selected values of the bias gap: (A) g = 0.7; (B) g = 0.9.

The evolution of the formation of inter-strand native contacts for the -sheet model is reflected in the snapshots in Fig. 7 at selected times that display typical folding events for: g = 0.9 (Fig. 7 A) and g = 1.1 (Fig. 7 B). As mentioned above, for g = 0.9 native contacts start forming at the ends of the chain as can be clearly seen in the structure at t* = 32. The native state is observed at t*  100. For g = 1.1, the chain collapse and the formation of native contacts are concurrent. An example of the collapsed chain is shown at t* = 32. A native-like structure is observed at t* = 63. A well-ordered native -sheet appears to exist at t*  100. 

View larger version (34K): [in this window] [in a new window]   FIGURE 7   Snapshots of the kinetic chain conformation of the -sheet at selected times, t* = 0.2, 32, 63, and 100 for the bias gaps: (A) g = 0.9; (B) g = 1.1.

Fig. 8 shows snapshots for the -clip at selected times as the chain evolves toward the native state. For g = 0.7 in Fig. 8 A, the four snapshots are taken after F^*_R reaches its minimum as seen in Fig. 2 B. The structure at t* = 100 corresponds to the collapsed chain at the minimum value of F^*_R. It exists for a very short time. This collapsed chain is a transient in protein folding (Silow and Oliveberg, 1997) and can be regarded as a kinetic disordered globule because it has a Q value, Q  0.6 (see Fig. 3 B), that is similar to that of a thermodynamic disordered globule (Jang et at., 2002) but with a smaller R. The existence of a disordered compact structure is also reported in the kinetic study of a model three-helix bundle protein (Zhou and Karplus, 1997b, 1999). At t* = 1.6 × 10³ and at t* = 1.3 × 10⁴, we see the structures of two typical intermediates; these are intermediate I₁ and intermediate I₂, respectively (see below). These intermediates are found in the strand extension process after the initial chain collapse and have very long lifetimes. An example of the native structure is shown at t* = 10⁵. For g = 0.9, the evolution of the chain toward the native state is totally different. No kinetic disordered globule state is observed when F^*_R goes to the minimum at t*  63. Instead, a slightly bent -sheet structure occurs as shown in Fig. 8 B at t* = 63. This -sheet-like structure is a transient in folding trajectories toward the -clip. Examples of the intermediate I and the native structure are observed at t* = 160 and t* = 10⁴, respectively.

View larger version (38K): [in this window] [in a new window]   FIGURE 8   Snapshots of the kinetic chain conformation of the -clip for two values of the bias gap: (A) g = 0.7 at times t* = 100, 1.6 × 103, 1.3 × 104, and 105; (B) g = 0.9 at times t* = 25, 63, 160, and 104.

Fig. 9 shows snapshots for the -twist at selected times as the chain evolves toward the native state. The behavior displayed in Fig. 9 A for g = 0.7 is similar to that shown in Fig. 8 A for the -clip. The collapsed chain at the minimum value of F^*_R is shown at t* = 63; this is a transient in protein folding (Silow and Oliveberg, 1997) and can be regarded as a kinetic disordered globule. The kinetic disordered globule has more non-native interactions than the thermodynamic disordered globule giving a smaller R. Due to the disordered nature of this state, no native-state secondary structure is observed. As time increases, the native contacts begin to form as part of the strand extension process, leading to the native structure. This process is very slow because the collapsed chain undergoes a subtle conformational change to the native state via intermediates that have long lifetimes. Examples of two typical intermediates, the intermediate I₁ and the intermediate I₂ are observed at t* = 250 and at t* = 1.3 × 10³, respectively. The native structure is observed at t* = 4 × 10⁴. For g = 0.9 in Fig. 9 B, no kinetic disordered globule state is observed when F^*_R goes to the minimum at t*  63. (This was also the case for the -clip at g = 0.9). Instead, a slightly bent -sheet structure is observed at t* = 63. The structures of the intermediate I (see below) are shown at t* = 320 and at t* = 3.2 × 10³. The native structure is observed at t* = 10⁴.

View larger version (42K): [in this window] [in a new window]   FIGURE 9   Snapshots of the kinetic chain conformation of the -twist for two values of the bias gap: (A) g = 0.7 at times t* = 63, 250, 1.3 × 103, and 4.0 × 104; (B) g = 0.9 at times t* = 63, 320, 3.2 × 103, and 104.

Intermediates and folding pathways

It is of interest to examine the intermediates that occur as the folding progresses in more detail. The contour graph in Fig. 10 A shows the population distribution of the fraction of native contacts, Q, as a function of reduced time on a logarithmic scale, log(1 + t*), for the -sheet at g = 0.9. The population distribution of Q is calculated at selected logarithmic time intervals of log(1 + t*) = 0.1 for all 100 trajectories. The graph clearly shows highly populated states as those enclosed by many contour lines. The -sheet at g = 0.9 is a fast-folding protein, and the folding trajectories clearly show a fast-track or bottleneck pathway of kinetic folding between two highly populated states: the initial random coil state at t* < 10 and Q  0.08 and the native state at t* > 100 and Q  1. More evidence of fast-track folding is shown in Fig. 10 B, which shows the population distribution of F^*_R versus Q averaged over all 100 trajectories and all times. The peaks indicate the highly populated states. The peak at F^*_R  0.8 and Q  0.1 corresponds to the initial random coil state and the other peak at F^*_R  0 and Q  1 corresponds to the native state. The kinetic folding mechanism for the -sheet at g  0.9 is represented in Fig. 10 C.

View larger version (29K): [in this window] [in a new window]   FIGURE 10   (A) The population distribution of the fraction of native contacts, Q, as a function of the reduced time in logarithmic scale, log(1 + t*); (B) the population distribution of the normalized squared radius of gyration, F*R, versus the fraction of native contacts, Q; (C) the kinetic folding mechanism for the -sheet at g = 0.9.

The contour graph for the distribution of Q versus log(1 + t*) for the -clip at g = 0.7 is shown in Fig. 11 A, and the corresponding population distribution of F^*_R versus Q from all 100 trajectories is shown in Fig. 11 B. In the contour graph of Fig. 11 A, the random coil state is observed in the range 0 < log(1 + t*) < 0.8 and 0 < Q < 0.1. Two long-lived populations are observed in the range 2.6 < log(1 + t*) < 3.4 at Q  0.76 and 3.4 < log(1 + t*) < 5 at Q  0.8; these correspond to the intermediates I₁ and I₂, respectively. Because the native state is located at t*  10⁵ as shown in Fig. 5 A for g = 0.7 and the results are shown within the time scale up to log(1 + t*) = 5, the population distribution of Q for the native state is not observed in the contour graph. In Fig. 11 B, the peak at F^*_R  1 and Q  0.76 corresponds to the intermediate I₁, which can be regarded as a partially-ordered globule as shown in the thermodynamic study of the same model (Jang et al., 2002). It has a small value of R, a large value of Q, and the rolled or kinked structure seen in Fig. 8 a at t* = 1.6 × 10³. This partially-ordered globule state occurs because the chain becomes trapped in a conformation dominated by the attraction between non-native contacts, thus disturbing the tendency to order into the native state. Note that the intermediate I₁ is highly populated at t*  10³ and has a value F^*_R  1, which is much smaller than the value of F^*_R = 0.4 at t*  10³ seen in Fig. 2 B, which is the average value of F^*_R for all trajectories. This suggests that not all 100 trajectories fold via the intermediate I₁. In fact, 40% of the trajectories are trapped in the intermediate I₁, and the rest of the trajectories fold via the other intermediate, I₂, which is located at F^*_R  0.2 and Q  0.8 in Fig. 11 B. The peak for intermediate I₂ is obscured by the peak of the highly populated native state at F^*_R  0 and Q  0.82. The intermediate I₂ is an ordered globule; it has a larger R than the intermediate I₁ and has almost the same Q as the native state. The intermediate I₂ has a loosely packed and slightly bent structure as seen in Fig. 8 A at t* = 1.3 × 10⁴. There are two possible folding pathways for the collapsed chain from the state in which F^*_R is at its minimum. The first is that it could further collapse to become trapped in the intermediate I₁, which can fold either directly into the native state or via the intermediate I₂ by elongating its strands. The other pathway is that it can fold via the intermediate I₂ to the native state. The kinetic folding scheme for the -clip at g = 0.7 is represented in Fig. 11 C.

View larger version (24K): [in this window] [in a new window]   FIGURE 11   (A) The population distribution of the fraction of native contacts, Q, as a function of the reduced time in logarithmic scale, log(1 + t*); (B) the population distribution of the normalized squared radius of gyration, F*R, versus the fraction of native contacts, Q; (C) the kinetic folding mechanism for the -clip at g = 0.7.

The folding kinetics for the -clip at g = 0.9 is very different than that at g = 0.7. Fig. 12 A shows the contour plot of the population distribution of Q versus log(1 + t*). The initial random coil state is located at 0 < log(1 + t*) < 0.6 and 0.05 < Q < 0.1. There is a short-lived intermediate I located at 0.75 < Q < 0.85 and 2.1 < log(1 + t*) < 2.7. The highly populated native state is located at t*  10⁴ with Q  0.9. The properties of intermediate I are the same as those of the intermediate I₂ for the g = 0.7 model. In Fig. 12 B, it is located at F^*_R  0.2 and Q  0.8, but its peak is obscured by the native peak at F^*_R  0 and Q  0.9. The intermediate I is an ordered globule with a native-like secondary structure as seen in Fig. 8 B at t* = 160. Almost 80% of the trajectories become trapped in the intermediate I, and the rest fold directly to the native state via a fast-track folding pathway. The kinetic folding scheme for the -clip at g = 0.9 is represented in Fig. 12 C.

View larger version (29K): [in this window] [in a new window]   FIGURE 12   (A) The population distribution of the fraction of native contacts, Q, as a function of the reduced time in logarithmic scale, log(1 + t*); (B) the population distribution of the normalized squared radius of gyration, F*R, versus the fraction of native contacts, Q; (C) the kinetic folding mechanism for the -clip at g = 0.9.

The folding kinetics for the -twist at g = 0.7 is qualitatively similar to that of the -clip at g = 0.7 (not shown in figures). The kinetic folding of the -twist involves two intermediates I₁ and I₂. However, unlike the -clip at g = 0.7 the intermediate I₁ is short-lived and relatively unpopulated in the F^*_R versus Q plane. Approximately, 20% of the trajectories are trapped in intermediate I₁, most of the trajectories fold into intermediate I₂, and none of the trajectories fold via the fast-track folding pathway.