{'fraction_non_alphanumeric': 0.03121655369247235, 'fraction_numerical': 0.02160631466286122, 'mean_word_length': 4.27511173841449, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Universality classes in folding times of proteins

Molecular dynamics simulations in simplified models allow one to study the scaling properties of folding times for many proteins together under a controlled setting. We consider three variants of the Go models with different contact potentials and demonstrate scaling described by power laws and no correlation with the relative contact order parameter. We demonstrate existence of at least three kinetic universality classes which are correlated with the types of structure: the alpha-, alpha--beta-, and beta- proteins have the scaling exponents of about 1.7, 2.5, and 3.2 respectively. The three classes merge into one when the contact range is truncated at a 'reasonable' value. We elucidate the role of the potential associated with the chirality of a protein.


INTRODUCTION
How do size and structure of a protein affect its folding kinetics is an interesting basic issue that has been debated in recent years. The size can be characterized by the number, N , of the amino acids that the protein is made of. The distribution of N across proteins stored in the data banks is peaked around N =100 (Cieplak & Hoang, 2000) and all proteins with a large N , like titin (N ≈ 30 000), consist of many domains. There must be then a mechanism that prevents globular proteins from reaching much larger sizes. We have argued (Cieplak & Hoang, 2000) that this is provided by the function of the protein which requires adoption of a specific conformation. Folding into it becomes increasingly difficult when N becomes larger and larger. The sizes of proteins are substantially smaller than those of the DNA molecules whose coding function does not depend on the shape. The native structure of a protein, on the other hand, is believed to be a decisive factor in its folding mechanism (Baker, 2000;Takada, 1999).
A simple parameter that is used to characterize the structure of the protein is the relative contact order, CO, (Plaxco et al., 1998) defined as average sequence distance between two aminoacids that interact with each other, i.e. form a contact, in the native state: where ∆ ij is 0 if the amino acids i and j do not form a contact and 1 otherwise. The relative contact order parameter is small for α-proteins in which all secondary structures consist of the α-helices because the hydrogen bonds in the helices correspond to |i − j| = 4. On the other hand, β-proteins tend to have larger CO because the β strands that form a sheet often involve amino acids which are quite distant along a sequence.
In their seminal 1998 paper (Plaxco et al., 1998) (paper I), Plaxco, Baker, and Simons have argued that folding rates correlate with CO but do not with N . Their argument was based on analysing experimental data on short proteins that were available in the literature. Their conclusion was reinforced in the 2000 paper (Plaxco et al., 2000) (paper II) by Plaxco, Simons, Ruczinski, and Baker in which the compilation of the kinetic data involved a larger set of proteins, including those that were considered in paper I. The later data were also restricted to a much narrower temperature range of between 20 and 25 o C. Their results for the folding times (i.e. the inverses of the folding rates) are represented in Figure 1 as a function of N (on the logarithmic scale). For the purpose of further discussion, we have divided the data into three classes: α-proteins, β-proteins, and α-β-proteins. The α-proteins are easily seen to be the fastest folders but clearly all of the data points are scattered all over the plane of the figure.
This random looking pattern of the data may, however, be only apparent since the plot might involve mixing distinct classes of proteins that perhaps should not be compared together. Figure 2 indeed hints at such a possibility as the splitting into the α-, β-, and α-β-structural classes reveals some patterns. These patterns are shown in different time windows -the α-proteins are in the window of much shorter times. There is a growing trend for the β-proteins and, if one disregards one outlayer, also for the α-β-proteins. The data for the α-proteins, however, are puzzling since if they do show an overall trend then it would be downwards, i.e. the bigger the N , the shorter the folding time which defies a common simple expectation to observe the opposite.
The combined data show a strong correlation with the CO parameter. When the data are split into the three structural classes, as shown on the right-hand side of Figure 2, then the correlation remains strong for the α-and α-β-proteins. However, in the crucial test case of the β-proteins (the right-bottom panel of Figure 2) four proteins have nearly the same CO and yet substantially different folding times. Thus there are some unsettling issues in our understanding of the experimental data that would be desirable to solve.
Theoretical modeling in simplified models, despite its well known general shortcomings, is expected to be a tool of help to identify possible trends because a unified approach can be applied to many different proteins. In this paper, we consider 51 proteins: 21 of the α-β kind with N between 29 and 162, 14 of the α-proteins with N between 35 and 154, and 16 β-proteins with N between 36 and 124. This set contains the 21 proteins, used in Figures 1 and 2, that were considered by Plaxco et al. All of the 51 proteins are modelled in three different ways and studied by the techniques of molecular dynamics. Even though the three models are all coarse grained and of the Go type (Abe & Go, 1981;Takada, 1999) they have very different kinetic and equilibrium properties when used for a particular protein. The variations between the models do lead to some differences in scaling properties of certain parameters, such as the temperature of the fastest folding, T min , or the thermodynamic stability temperature, T f , but they all agree on a power law dependence of the folding time, and on the lack of any correlation of t f old with CO. The problems with the experimental results on the N dependence may be related to the lack of the temperature optimization. The folding time often depends on the temperature, T , and making choices on the temperature to study kinetics may affect the outcome of the measurement. We argue that a demonstrable trend might arise when all data are collected at T min which needs to be determined for each protein individually.
The theoretically derived lack of correlation of t f old with CO seems to be a more difficult issue. One may just dismiss it as characterizing not real life but an approximate model. On the other hand, the essence of the Go models is that they are based on the native topology. Thus if such geometry sensitive models do not 'care' about the contact order then what models would? We leave it as an open question and this paper may be just considered to be a report on what are the properties of three different Go-like models. Notice, however, that once a Go model is constructed its contacts are well defined and the kinetics are studied in the context of such a definition whereas assignment of contacts in experimental systems is subjective. It should be pointed out that the contact order in the Go models is actually quite important but not for the overall folding time -it is the primary factor that governs the succession of events during folding (Unger, 1996;Hoang & Cieplak, 2000a;Hoang & Cieplak, 2000b;Cieplak et al., 2002a;Cieplak et al., 2002b;Erman, 2001). In other words, what is important for folding of a protein is the full "spectrum" of the relevant values of the sequence distances |i−j| and not just their average value. A similar point has been argued within a host of models in references (Galzitskaya and Finkelstein, 1999;Alm and Baker, 1999;Munoz and Eaton, 1999;Du et al., 1999;and Plotkin and Onuchic, 2000).
The power law dependence described by eq. 2 has been proposed by Thirumalai (Thirumalai, 1995) and then demonstrated explicitly for several types of lattice models (Gutin et al., 1996;Zhdanov, 1998;Cieplak et al., 1999). On the other hand, a number of theories and a recent simulation of 18 proteins (away from the optimal folding condition) by Koga and Takada (Koga & Takada, 2001) suggest a power law dependence for barrier heights on N and hence an exponential dependence of t f old on N (Takada & Wolynes, 1997;Finkelstein & Badredtinov, 1997;Wolynes, 1997). Thus the issue of scaling remains unsettled not only experimentally but also theoretically.
Recently (Cieplak & Hoang, 2001) we have demonstrated the power law dependence for one variant (of the three studied here) Go-like models when applied to 21 proteins which were mostly of the α-β kind. The resulting exponent λ turned out to be equal to 2.5 ± 0.2. In this particular variant of the Go model, the native contact interactions were restricted to a cut-off value of 7.5Å and the contact potential was described by the Lennard-Jones form.
Here, we extend such studies to the other two kinds of proteins, α and β, and arrive at a similar value of λ. However, when the model is made significantly more realistic by considering the range of the native contact interactions as a variable quantity, then we arrive at a richer picture. We show that the three classes of tertiary structures also correspond to three different kinetic universality classes. The α-proteins come with λ of around 1.7 (the result obtained previously (Cieplak & Hoang, 2001) for decoy helical structures), the β-proteins are characterized by λ close to 3.2, and the α-β-proteins have λ near 2.5. These values do not depend on whether the contact potential are Lennard-Jones or of the 10-12 form so they are truly a reflection of the native topology. The power law trends are pretty evident when the folding times are determined at T min but harder to see otherwise. In these studies, the range of the contact interactions has been determined based on the van der Waals radii of the atoms (Tsai et al., 1999). Another realistic item that we implement is the chirality potential -a term which is responsible for folding to a conformation of the correct native chirality. This term affects the kinetics but we show it not to affect values of the exponent λ.
The growth of t f old with N indicates increasingly deteriorating folding conditions. Our studies of scaling of T min and T f indicate that asymptotically T f becomes substantially lower than T min which signifies an onset of slow glassy kinetics before the system is near the native conformation. This adds to the deterioration of foldability and suggests the limitation in the observed values of N . The three models considered here have T min and T f varying as a function of N in different ways, though they agree asymptotically. Among the three models, the Lennard-Jones contact potential with the variable R c appears to have the most appealing kinetic properties in that it leads to a very good foldability for a small N . This should be our simple model of choice in future studies. However, the issue of the scaling trends needs now to be studied in models that reach beyond the Go approximation and in experiments with a protocol that involves optimization.

A. The Hamiltonian
An input for the construction of the Go model is a PDB file (Bernstein et al., 1977) with the coordinates of all atoms in the native conformation. The coordinates are used to determine the length related parameters of the model. Whereas all energy and temperature related parameters are expressed in terms of a common unit -ǫ. We model 51 proteins. In addition to the proteins listed in the caption of Figure 1, we also consider 1cti ( We consider several variants of the Go models. In each case, the Hamiltonian consists of the kinetic energy and of the potential energy, E p ({r i }), which is given by The first term, V BB is the harmonic potential which tethers consecutive beads at the equilibrium bond length, d 0 , of 3.8Å. Here, r i,i+1 = |r i − r i+1 | is the distance between the consecutive beads and k = 100ǫ/Å 2 , where ǫ is the characteristic energy parameter corresponding to a native contact.
The native contacts are defined either through the distances between the C α atoms or through an all-tom consideration. The first choice, used by us previously (Hoang & Cieplak, 2000a;Hoang & Cieplak, 2000b;Cieplak & Hoang, 2001), is to take a uniform cut-off distance, R c , of 7.5Å, below which a contact is said to be present. In the second choice, used here in most cases, all the heavy atoms present in the PDB file are taken into account. Specifically, a pair of aminoacids is considered to form a contact if any pair of their non-hydrogen atoms have a native separation which is smaller than 1.244 (R i + R j ), where R i are the van der Waals radii of atom i, as listed in ref. (Tsai et al., 1999). This critical separation corresponds to the point of inflection of the Lennard-Jones potential. Figure 3 shows the distribution of the effective contact ranges as obtained for an N =162 protein T4 lysozyme with the PDB code 102l which consists of 10 α-helices and 3 β-strands. There are 339 native contacts in this case and they range in value between 4.36 and 12.80Å. It is clear that truncating this distribution at whatever "reasonable" value, which is often taken to be in the range between 6.5 and 8.5Å would result in a substantial removal of the relevant interactions. Thus insisting on a uniform cutoff value is expected to have noticeable dynamical effect.
We consider two variants of the interactions in the native contacts. The first variant is the 6-12 Lennard-Jones potential where the sum is taken over all native contacts. The parameters σ ij are chosen so that each contact in the native structure is stabilized at the minimum of the potential, and σ ≡ 5Å is a typical value. The second variant is the 10-12 potential where r (n) ij coincides with the native distance. This potential is frequently used to describe hydrogen bonds (Clementi et al., 2000). For each pair of interacting amino acids, the two potentials have a minimum energy of −ǫ and are cut off at 20Å. The non-native interactions, V N ON , are purely repulsive and are necessary to reduce the effects of entanglements. They are taken as the repulsive part of the Lennard-Jones potential that corresponds to the minimum occurring at 5Å. This potential is truncated at the minimum and shifted upward so that it reaches zero energy at the point of truncation.
The final term in the Hamiltonian takes into account the chirality. Natural proteins have right handed helices but a Go model as described above involves chiral frustration: one end of a helix may want to fold into a right handed helix and another into a left handed one and "convincing" one end to agree with the twist of the other takes time and delays folding. Such a frustration would not arise naturally. In order to prevent it, we add a term which favors the native sense of the overall chirality at each location along the backbone. A chirality of residue i is defined as where v i = r i+1 − r i . A positive C i corresponds to righthanded chirality. Otherwise the chirality is left-handed.
The values of C i are essentially between −1 and +1. The distribution of C i in 21 α-β-proteins considered in this study is shown in Figure 4. It is seen to be bimodal. The values in the higher peak correspond to locations within the helical secondary structures. The chiral part of the Hamiltonian is then given phenomenologically by where Θ is the step function (1 for positive arguments and zero otherwise), C N AT i is the chirality of residue i in the native conformation, and κ is taken, in most cases, to be equal to ǫ. However, a criterion for selection of its proper value remains to be elucidated. The idea behind this particular form of V CHIR is that when the local chirality agrees with the native chirality then there is no effect on the energy. On the other hand, a disagreement in the chirality is is punished by a cost which is quadratic in chirality.
V CHIR has the strongest effect on the helical structures. However, it affects the sense of a twist of the whole tertiary structure. The chirality term enhances the dynamical bias towards the native structure during the folding process and helps avoiding non-physical conformations such as left-handed helices. V CHIR is a four-body potential. In this respect this term is similar to potentials that involve dihedral angles (Veitshans et al., 1997;Clementi et al., 2000;Settanni et al., 2002). The dihedral terms enhance stability of a model of the protein but usually have no bearing on the chirality (Veitshans et al., 1997) unless they involve directly the values of native dihedral angles (Clementi et al., 2000;Settanni et al., 2002).

B. The time evolution
The time evolution of unfolded conformations to the native state is simulated through the methods of molecular dynamics as described in details in (Hoang & Cieplak, 2000a;Hoang & Cieplak, 2000b) (see also (Cieplak et al., 2002a;Cieplak et al., 2002b) in the context of the Lennard-Jones contact potentials. The beads representing the amino acids are coupled to Langevin noise and damping terms to mimic the effect of the surrounding solvent and provide thermostating at a temperature T . The equations of motion for each bead are where m is the mass of the amino acids represented by each bead. A similar approach in the context of proteins has also been adopted in references (Guo and Thirumalai, 1996;Berriz et al., 1997;and Eastman and Doniach 1998). The specificity of masses has turned out to be irrelevant for kinetics (Cieplak et al., 2002a) and it is sufficient to consider masses that are uniform and equal to the average amino acidic mass. F c is the net force due to the molecular potentials and external forces, γ is the damping constant, and Γ is a Gaussian noise term with dispersion √ 2γk B T . For both kinds of the contact potentials, time is measured in units of τ ≡ mσ 2 /ǫ, where σ is 5Å. This corresponds to the characteristic period of undamped oscillations at the bottom of a typical 6-12 potential. For the average amino acidic mass and ǫ of order 4kcal/mol, τ is of order 3ps. According to Veitshans et al. (Veitshans et al., 1997), realistic estimates of damping by the solution correspond to a value of γ near 50 m/τ . However, the folding times have been found to depend on γ in a simple linear fashion for γ > m/τ (Hoang & Cieplak, 2000a;Hoang & Cieplak, 2000b;. Thus in order to accelerate the simulations, we work with γ = 2m/τ but more realistic time scales are obtained when the folding times are multiplied by 25. The equations of motion are solved by means of the fifth order Gear predictor-corrector algorithm (Gear, 1971) with a time step of 0.005τ .
The magnitude of the viscous effects, as controlled by the parameter γ, has to be sufficiently large so that the scenarios of the folding events are not dominated by the inertial effects. Otherwise the scenarios would depend on the spacial and not on the sequencial separation between the amino acids. Figure 5, for crambin as an illustration, shows that even though our value of γ of 2 is reduced compared to the values that are expected to be realistic it already corresponds to sufficiently strong damping with the minimal inertial effects. Figure 5 gives average first times needed to establish contacts separated by the sequence length |i − j| for three values of γ: 2, 12, and 24 m/τ . To the leading order, the times to establish the contacts (and also the folding times) are linear functions of γ so one can show them together by proper rescaling. Furthermore, the whole pattern of the events is insensitive to the value of γ. Starting with this figure, we adopt the convention that the symbol sizes give measures of the error bars in the quantity that is plotted.
The folding time is calculated as the median first passage time, i.e. the time needed to arrive in the native conformation from an unfolded conformation. It is estimated based on between 101 and 201 trajectories. T min is defined as a temperature at which t f old has a minimum value when plotted vs. T . For small values N , the U-shaped dependence of t f old on N may be very broad and then T min is defined as the position of the center of the U-shaped curve. The simplified criterion for an arrival in the native conformation to be declared is based on a simplified approach in which a protein is considered folded if all beads that form a native contact are within the cutoff distance of 1.5σ ij or 1.2r (n) ij for the 6-12 and 10-12 potentials respectively.
The stability temperature T f is determined through the nearly equilibrium calculation of the probability that the protein has all of its native contacts established. T f is the temperature at which this probability crosses 1 2 . The calculation is based on least 5 long trajectories that start in the native state in order to make sure that the system is in the right region of the conformation space. It should be noted that, in the literature, the frequently used estimate of the folding temperature is determined through the position of the maximum in the specific heat. This yields a T ′ f which is typically larger than T f . Our probabilistic interpretation has the disadvantage of being dependent on the precise definition of what constitutes the native basin (and thus only the approximate location of T f is of relevance) but it has the advantage of relating only to the native basin and not to any other valleys in the phase space. In most of our systems, T f is found to be comparable to T min , while both of them are always lower than T ′ f . Furthermore, in most cases, even though when T min is found to be higher than T f , the folding times at T f are comparable to those at T min which indicates that the model is unfrustrated in the conventional sense. Only in some very few cases, the folding times at T f are excessively long to be determined in our simulation. This behaviour probably corresponds to a structural frustration (Clementi et al., 2000) embedded in the native conformation.
An alternative to the contact-based criterion for folding is to provide a more precise delineation of the native basin as in ref. (Hoang & Cieplak, 2002b) or relate the criterion to a cutoff in the value of the RMSD distance away from the native conformation. These approaches are illustrated in Figure 6 which shows the dependence of the folding time, t f old , vs. T for a synthetic α-helix (H16 of reference (Hoang & Cieplak, 2002a)) and β-hairpin (B16 of the same reference) that both consist of 16 monomers. Whichever criterion for folding is used, the folding curves are U-shaped and the non-zero chirality term extends the region of the fastest folding both towards the low and high temperature ends. For the hairpin, the effect is smaller but still clearly present.
When it comes to model proteins, we used only the contact-based folding criterion. An illustration of the role of the chirality potential is provided in Figure 7 for crambin (N =46, the PDB code 1crn) which is a protein of the α-β type. The top panel, for R c = 7.5Å, shows that the shortest time of folding is somewhat reduced by V CHIR but the biggest impact is on the range of temperatures at which folding is optimal, almost by the factor of 2, especially in the low T regime. For the β proteins, the effect of the chirality potential is generally smaller. For the SH3 domain coded 1efn the change due to V CHIR is hard to detect (not shown) but for the I27 globular domain of titin, coded 1tit, it is quite substantial on the low T side of the curve (Figure 8). We conclude that incorporation of the chirality term in the Hamiltonian appears to reduce structural frustration in these models and thus makes the models more realistic. For all of the results presented here from now on (except for Figure 12), the chirality term is included.
Another simple way to enhance the realism of the Go models is suggested by Figure 3: calculate the range of the contact potential instead of taking one uniform cutoff value. When we compare the case of the Lennard-Jones contact potential with the uniform or variable R c then the nature of the effect on the kinetics strongly depends on the protein. For instance, for the protein 1crn ( Figure  7, bottom panel) there is essentially no difference. On the other hand, a dramatic narrowing of the U-curve is observed for 1tit (Figure 8).
On switching the 6-12 potential to the 10-12 potential all of the kinetic U-curves become substantially narrower (Figures 7 and 8). This is related to the fact that the potential well corresponding to the 10-12 potential is narrower which makes folding a task that requires more precision. Note, that the two potentials have the same energy (−ǫ) at the minimum so the temperature scale are comparable.
We have demonstrated that there are many ways to construct variants of the Go models and they all come with distinctive folding characteristics.

RESULTS
A. The 6-12 potential with the variable contact range Figure 9 shows the median values of t f old at T min for the Lennard-Jones contact potential when the presence of the native contact is determined through the van der Waals sizes of the atoms (and with the chirality term included). Figure 9 data divides the data into the three structural classes. There are a few outlayers (one is the 1aps protein which appears to be a poor folder also experimentally) but basically there are clear linear trends on the log-log scale which indicates validity of the power law, eq. (2). The values of the exponents 1.7 for the α-proteins and 3.2 for the β-proteins agree with those found for decoy structures (Cieplak & Hoang, 2001). The decoy structures were constructed from homopolymers and the contact range was not variable due to the lack of atomic features in the decoys. Figure 10 replots the same data together to indicate that the trends identified in the classes are identifiably distinct. Thus the structural classes also correspond to the kinetic universality classes. Figure 11 shows data equivalent to those on Figure 9 but now the folding times are determined at T f , as an example of a situation that may be encountered away from the optimal conditions. The data points show a much larger scatter away form the trend identified at T min . The optimal trend seems still dominant but it is so much harder to see. This should be analogous to results obtained experimentally.
It is interesting to figure out what is the effect of the chirality potential on the scaling results. Figure 12 refers to the α-proteins and it compares the case of κ = 0 to κ = ǫ. Proteins with small values of N are not sensitive to the value of κ but for N >≈ 50 taking the chirality into account accelerates the kinetics quite noticeably. The 'asymptotic' scaling behavior remains unchangedthe exponent λ of 1.7 is valid for both cases, though a somewhat larger value for κ=0 cannot be ruled out (but certainly not as large as 2.5). We have checked that the data points for κ = 2ǫ, though corresponding to a bit faster times than for κ = ǫ, are in practice indistinguishable from the latter in the scale of the figure. This observation suggests a behavior which saturates with a growing κ.
As pointed out in Ref. (Cieplak et al., 1999), the dependence of T f and T min on N may offer additional clues about the foldability at large N . Figure 13 suggests that the α-and α-β-proteins are excellent folders for small values of N since then T min is less than T f . T f appears to have no systematic trend with N but the data for T min suggest a weak growth, approximately proportional to log(N ). Around N of 50 the trend associated with T min crosses the average value of T f and from now on T f is lower than T min . This suggests that asymptotically the energy landscape of the system would be too glassy-like to sustain viable folding. Thus accomplishing folding would require breaking into independently folding domains domain or receiving an external assistance, e.g. from chaperons whereas our studies are concerned with individual proteins. Figure 13 also suggests that the β proteins behave somewhat differently since they exhibit no trend in T min in the range studied and already for small values of N T min exceeds T f . Nevertheless the differences between the three structural classes are minor because they all show a border line behavior: the proteins in the range up to N =162 are not excellent but just adequate folders, at least in this model.
It is interesting to point out that neither t f old nor the characteristic temperatures indicate any demonstrable correlation with the relative contact order defined in eq. 1. This is shown in Figure 14: for a given value of CO we find systems both with long and short folding times or both high and low values of T min . B. The 10-12 potential with the variable contact range We now check the stability of our results against the change in the form of the contact potential with the same characteristic energy scale. Figure 15 shows that when the Lennard-Jones potential is replaced by the 10-12 potential, with keeping all other Hamiltonian parameters intact, the scaling trends for T f old are consistent with those displayed in Figure 9 and confirm the existence of the three universality classes. Figure 16 suggests that the 10-12 systems are also border line in terms of the positioning of T min vs T f but the weak growing trends for the α-and α-β-proteins are gone. The lack of correlations with the relative contact order also holds for the 10-12 potential (not shown).
C. The 6-12 potential with R c = 7.5Å We now return to the Lennard-Jones potential and make the drastic, as evidenced by Figure 3, change that only those native contacts are considered whose range does not exceed 7.5Å. The resulting data are shown in Figure  17. The top panel indicates that λ of about 2.5 is still consistent with the trend obtained. However, λ of 1.7 is quite off the mark for the α-proteins. The exponent of 3.2 for the β-proteins is not ruled out but the scatter in the data points is bigger than in the bottom panel of Figure  9. Taken together with the results for the α-proteins, the most likely conclusion is that the fixed, and invasive, cut off in the contact range looses the ability to distinguish between the structural classes and all such models of the proteins would be characterized by a single exponent λ of 2.5 as found in ref. (Cieplak & Hoang, 2001). This is illustrated in Figure 18 where the data corresponding to various structural classes are displayed together. They seem to be consistent with just one trend. Figure 19 shows T min and T f for the case with R c = 7.5Å. It suggests that among the three models studied here, the one with the cut off in the contact range is the worst kinetically because the gap between the band of values of T min and the band of values of T f is the largest. This indicates that precise values of the contact range are important in the task of putting pieces of a protein together in the folding process. Also in this model, there is no correlation with the relative contact order parameter.

DISCUSSION
We have studied 3 variants of the Go model through the molecular dynamics simulations and demonstrated the power law dependence of the folding time on N and lack of dependence on CO. Furthermore, the models with the variable contact range allow one to identify (at least) three kinetic universality classes corresponding to three different values of the exponent λ. The lowest exponent found for the α-structures is consistent with the widely held belief that the α-helices are structures that are optimal kinetically (Micheletti et al., 1999;Maritan et al., 2000). The scaling behavior of T min and T f , taken together with the increasing t f old suggests an asymptotic emergence of a glassy behavior. As a technical improvement, we have highlighted benefits of introducing the chirality potential.
Recently, Koga and Takada (Koga & Takada, 2001) have also studied scaling of t f old in proteins approximated by the Go model. They have considered the 10-12 potential that was augmented by potentials which involved the dihedral angles (but no chirality). They have determined the folding temperature through the maximum in the specific heat. Their studies at T ′ f , done for 18 proteins with N in the range between 53 and 153, suggest a t f old that exponentially depends on the relative contact order multiplied by N 0.6 . It is thus interesting to check on this conclusion in the framework of our approach. Figure 20 shows log(t f old ) vs. CO×N 0.6 for our best model, i.e. for the Lennard-Jones contact potential with variable contact range. It is clear that the data at T min (the left panels) show significantly less scatter than at T f (the right panels) so the distinction between the power law and the exponential function is certainly not due to considering different temperatures. Figure 20 does suggest a correlation with CO×N 0.6 (the data plotted vs. N 0.6 without the CO factor have a similar appearance indicating the irrelevance of CO in such theoretical studies) and Koga and Takada quote a correlation level of 84% for their data. It is not very easy to distinguish between the power law and the exponential dependencies without a significant broadening of the range in the values of N . Figure 21 shows the data of Figure 9 redisplayed on the log -linear scale. The exponential trends, t f old ∼ exp(N/ξ), cannot be ruled out and the correlation levels are 75%, 94%, and 95% for the α-β, α, and β structural classes respectively whereas the corresponding values for the log-log plots are 81%, 97%, and 94%. Even though the power law fits appear better (or, in the case of the β proteins about the same) the important point is that the exponential fits also suggest existence of the three different kinetic universality classes since the characteristic values of the ξ parameter, as displayed in the Figure, are clearly distinct. Our trends displayed in Figure 9 seem much less scattered than those shown in Figure 20, especially in the right hand panels of Figure 20. However, while we argue in favour of the three universality classes and then the power laws, we see a need for further studies and better understanding of these issues.
It has been found recently (Cieplak & Hoang, 2002c) that the kinetics of Go models are very sensitive to the selection of what constitutes the proper set of the native contacts. For instance, if one declares a uniform cutoff range, R c , between the C α atoms for making a contact, then the dependence of t f old on R c is strong and non-monotonic. Koga and Takada declare the contact as occurring if two non-hydrogen atoms in a pair of amino acids are in a distance of less than either 5.5Å or 6.5Å (and it is stated that the results are stable with respect to this choice). Our definition of the contacts, on the other hand, involves the atomic sizes which yields a different contact map and leads to different folding times.
The basic unsolved question is why do the folding times in various Go models do not depend on the contact order even though the primary ingredient of any Go model is the geometry of the native state of a protein. One technical problem with the contact order is that the very notion of a contact is fairly subjective. Consider, for instance, the G protein -the PDB code is 1gb1 for the structure determined by NMR and 1pga for the crystallographic structure. When we make use of the van der Waals radii then we get CO = 0.239 for 1gb1 and 0.250 for 1pga. The alternative procedure is to consider two residues contacting if they contain non-hydrogen atoms within a distance of d. For d equal to 3, 4, 5, 6, 7 and 8 A, our procedure yields CO of 0.194, 0.220, 0.235, 0.252, 0.277, and 0.295 respectively (for the 1pga structure it is 0.257 if the cutoff of 6Å is used -i.e. not very different). Plaxco et al. (Plaxco et al., 1998;Plaxco et al., 2000) used the value of d = 6Å, and they quoted CO of 0.173 for this case. The notable difference from our value arises from the fact that in their calculation (Plaxco -private communication), all of the contacts made by the atoms (i.e. up to dozens for a pair of amino acids) contribute to the value of CO if the corresponding distance does not exceed d. Furthermore, the 'contacts' between consecutive residues (i.e. between i and i + 1) are taken into account. In our calculation, the shortest local contacts are of the i, i + 2 type. Note that the values of CO vary with d quite substantially (on the scale of the figures involved) and the value obtained at d = 6Å is about 45% larger than that quoted by Plaxco et al. The important point, however, is not that much what is the absolute value of CO but whether its correlation with the folding rate is sensitive to the choice of a specific definition of CO that is adopted. We have found that, quite remarkably, this correlation in the set of the experimentally studied proteins remains strong even when our procedure for the calculation of CO is used. We find that even though the scatter away from the trend is noticeably larger than when using the CO P -the values of CO quoted by Plaxco et al. -the correlations with CO remain robust and some dependence on CO develops in the case of the β-proteins. It is hoped that further interactions and iterations between theory and experiment will make the issues of size and contact order dependence more definitive. The notion of universality classes in proteins should play an important role in this process. Figure 1. Experimentally determined folding times based on tables compiled by Plaxco et al. (Plaxco et al., 2000). The solid circles, open hexagons, and stars are for the α-β-, α-, and β-proteins respectively. Figure 2. Experimentally determined folding times as split into three structural classes. The panels on the left hand side show the dependence on N whereas the panels on the right hand side show the dependence on the relative contact order parameter. Note that the time window in the middle panels are shifted by two orders of magnitude compared to the other panels. The top panels corresponds to the following α-β-proteins: 1div(56), 1gb1(56), 2ptl(63), 2ci2(65), 1aye(71), 1ubq(76), 1hdn(85), 2u1a(88), 1aps(98), and 2vik (126), where the number in brackets indicates the value of N (in the case of 2ci2 there are 19 more amino acids but their structure is undetermined). The middle panels correspond to the following α-proteins: 2pdd(43), 2abd(86), 1imq(86), 1lmb(92), 1hrc(104), 256b(106), and 1f63 (154). The bottom panels correspond to the following β-proteins: 1efn(57), 1csp(67), 1ten(89), and 1tit(89). The papers by Plaxco et al. (Plaxco et al., 1998;Plaxco et al., 2000) also contain data on several other proteins that are not shown here -we restrict ourselves only to the proteins that we study through simulations. (We had difficulties with the identification of the proper structure files for the remaining proteins). The subscript P in CO P signifies the criterion of Plaxco et al. (Plaxco et al., 1998;Plaxco et al., 2000) for a formation of a contact: two residues are considered to be contacting if they contain non-hydrogen atoms within the distance of 6Å. The symbols are as described in the caption of Figure 1. Figure 3. The distribution of the effective contact lengths in T4 lysozyme as determined by the procedure which is based on the van der Waals radii of the atoms. The shaded region corresponds to the contacts that would not be included if the cutoff of 7.5Å was adopted. Figure 4. The distribution of the chirality parameter C in 21 α-β-proteins studied. Figure 5. Times to establish contacts of a given sequence separation, |i − j| for crambin and for the indicated values of the damping constant γ. The times are rescaled so that k is equal to 1, 6, and 12 for γ equal to 2, 12 and 24 m/τ respectively and shown top to bottom. The symbols corresponding to γ = 12m/τ are reduced in size for clarity. The magnitude of the remaining symbols indicates the size of the error bars. The model used here corresponds to the Lennard-Jones contacts and the contacts are determined based on the van der Walls radii. The criterion for establishing a contact (for the first time) is based on whether the two beads come within a distance of 1.5σ ij of each other. This figure illustrates existence of second order effects in the dependence on γ because the rescaling by k brings the data points for a given event together but there is no strict overlapping. Figure 6. The dependence of the folding time on temperature for "synthetic" secondary structures of 16 monomers. The top two panels are for the α-helix system H16 and the bottom panel is for the β-hairpin B16. The dotted lines correspond to the chirality potential (with κ=1) included and the solid lines are for the case when it is not. In all cases, R c = 7.5Å. The top panel corresponds to the contact based criterion whereas the other panels is for the criterion based on the cutoff RMSD of 0.2Å. On the bottom panel, both curves include the chirality potential. Here, the solid (dashed) line is for the 6-12 (10-12) contact potential. The arrows indicate values of the folding temperature T f . The heavier (lighter) arrow is for the 6-12 (10-12) potential. Figure 8. The dependence of the folding time on temperature for models of the protein 1tit. The symbols are as in Figure 6: the thin solid line and the triangular data points are for R c = 7.5Å and no chirality; the dotted line with the square data points are for R c = 7.5Å and with the chirality; the thick solid line with the solid circular data points are for the locally calculated R c and the Lennard-Jones contact potential with the chirality; the dashed line with the open circular data points are for the similar case with the 10-12 potential. The arrows indicate the values of T f for the contacts of variable range: thick for the Lennard-Jones case and thin for the 10-12 case. Figure 9. The scaling of t f old with N for the 51 proteins as modeled by the 6-12 contact potential with the variable contact range. The data are split into the α-β-, α-, and β-proteins as indicated. The lines indicate the power law behavior with the λ exponent displayed in the right hand corner of each panel. The error bars in the exponent are of order ±0.2. The folding times are calculated at T min . The correlation levels of the points shown are 81%, 97% and 94% for the top, middle and bottom panels respectively. broken lines correspond to the slopes of 3.2, 2.5, and 1.7 respectively. The correlation level is 87%. Figure 11. Same as Figure 8 but the folding times are determined at T f instead at T min . The data points represented by the arrows indicate values which are significantly off the frame of the figure (for which only the lower bound of 30000 τ is known). The correlation levels are 83%, 88% and 77% for the top to bottom panels respectively. Figure 12. The role of the chirality potential on the folding times for the α-proteins. The hexagons are the data points shown in the middle panel of Figure 8 whereas the crosses correspond to the results obtained for κ=0. Figure 13. The values of T min and T f shown vs. N for the Lennard-Jones potential with the variable contact range. The data points are divided into the three structural classes. Figure 14. The dependence of t f old , T min , and T f on the relative contact order parameter for the Lennard-Jones contact potential with the variable contact range. The data symbols indicate the structural classes and are identical to those used in Figures 8, 9, 10, and 11. Figure 15. Same as in Figure 8 but for the 10-12 contact potential. The correlation levels are 88%, 98% and 91% from top to bottom. Figure 16. Same as in Figure 12 but for the 10-12 contact potential. Figure 17. Same as in Figure 8 but for the Lennard-Jones potential with R c = 7.5Å. The correlation levels are 83%, 91% and 93% for the top to bottom panels respectively. Figure 18. Same as in Figure 9 but for the cutoff of 7.5Å in the range of the contact potential. The solid line has a slope of 2.5. The correlation level for all of the points is 88%. Figure 19. Same as in Figure 12 but for R c = 7.5Å.   figure 9 redisplayed on the loglinear plane. The dashed lines indicate fits to the exponential law t f old ∼ exp(b/ξ) with the values of ξ shown in the right hand corner of each panel. The correlation levels are 75%, 94% and 95% for the top to the bottom panels respectively. The overall correlation level is 82% whereas for the power law fit it is 86%. The corresponding numbers for the 10-12 potential and the Lennard-