INTRODUCTION

TOP ABSTRACT INTRODUCTION NATIVE-STATE ANALYSIS NON-NATIVE INTERACTIONS VALIDITY OF G MODELS METHODS REFERENCES

Protein folding is one of the essential reactions in living systems. Recently, attention has focused on this reaction not only because of its fundamental role (Fersht, 1999), but also because of the interest in protein folding generated by the availability of many protein sequences from a rapidly increasing number of genomes and the realization that misfolded proteins are involved in disease (Dobson, 1999, 2001). Considerable progress has been made in achieving an understanding of the folding reaction by the use of simplified models (Bryngelson et al., 1995; Chan and Dill, 1998; Dinner et al., 2000). In particular, the problem posed by the "Levinthal Paradox" (namely that a polypeptide chain can find its unique native structure in spite of the very large number of possible denatured conformations) has been solved. It has been shown that a reasonable energy bias toward the native state can reduce the search of conformation space sufficiently for folding to take place on the experimental time scale (Karplus, 1997). One model, referred to as the "G model" or "G-type model" (Taketomi et al., 1975; Takada, 1999), has been widely used in studies of protein folding (Zhou and Karplus, 1999; Alm and Baker, 1999; Muñoz and Eaton, 1999; Galzitskaya and Finkelstein, 1999; Ozkan et al., 2001; Vendruscolo et al., 2001; Shimada et al., 2001). It is characterized by an energy function that replaces the nonbonded interactions (van der Waals and electrostatic terms) by attractive native-state contact energies; in some cases, non-native repulsions are also present (Zhou and Karplus, 1999; Shimada et al., 2001). Applications of G-type models include one-dimensional residue-based phenomenological descriptions of the folding reaction (Alm and Baker, 1999; Muñoz and Eaton, 1999; Galzitskaya and Finkelstein, 1999), lattice model calculations (Ozkan et al. 2001), and three-dimensional C or all-atom folding simulations by molecular dynamics (Zhou and Karplus, 1999) and Monte Carlo methods (Shimada et al. 2001); in the latter an excluded volume term is added to prevent collapse of the structure. For the phenomenological descriptions (Alm and Baker, 1999; Muñoz and Eaton, 1999; Galzitskaya and Finkelstein, 1999), the G-type model provides an essential simplification, which makes possible the replacement of the three-dimensional structure of the protein by a one-dimensional construct. For simulations in three dimensions (lattice and off-lattice), G-type potentials have the important property that the native state is a deep minimum, and that the potential surface corresponding to the non-native configurations is relatively smooth. There results a nearly ideal folding "funnel" leading to the native state, in contrast to the much rougher energy surface obtained with a more realistic molecular mechanics potentials (Duan and Kollman, 1998). As a consequence, trajectories calculated with G-type potentials take only on the order of nanoseconds to fold, instead of the experimental timescale which is microseconds or longer. This has made it possible to obtain statistically meaningful results for generic G-type models of proteins and polypeptide chains, and for models with native structures corresponding to those of specific proteins (Zhou and Karplus, 1999; Vendruscolo et al., 2001; Shimada et al., 2001).

Given the widespread use of G-type models, it is surprising that no direct tests have been made to determine whether they provide an accurate description of the protein energy surface. In this paper, we make such a test by comparing G-type model results with those obtained from an extensively validated effective energy function (EEF1) of the molecular-mechanics type, which combines a standard representation of the nonbonded van der Waals and electrostatic contributions to the energy with an implicit treatment of solvation (Lazaridis and Karplus, 1999). In an earlier paper (Paci et al., 2002b) we showed that the contact approximation with a cut-off radius of 5.5 Å for the EEF1 potential gives an excellent approximation to the total energy for native and non-native configurations of proteins, even when non-native interactions contribute significantly. It was also shown there that the inclusion of solven shielding is essential for the validity of the contact model. Here, we examine the validity of the most widely used form of the G-type model, which assumes a simple relation between the residue-residue interaction energies and the number of contacts. Results for the original G model (Taketomi et al. 1975), which uses a single parameter to relate a residue-residue contact to its interaction energy, are similar.

EEF1 and G-type models correspond to potentials of mean force; i.e., they represent the effective energy of the protein-solvent system for a given configuration of the protein in the presence of a canonically averaged solvent (Karplus and Shakhnovich, 1992). The EEF1 energy function can be decomposed into a sum of pairwise residue-residue interactions, E_IJ, where I and J correspond to the residues. For each geometry or for an ensemble of geometries, such as those representing the unfolding state (see Methods), we can write the EEF1 energy in the form

(1)

where we excluded N near-neighbor residue interactions (we use N = 2 in accord with several implementations of G-type models (Muñoz and Eaton, 1999; Vendruscolo et al., 2001; Shimada et al., 2001)); for N  2, the bonded terms make no contribution, and E_IJ can be written (see Methods)

(2)

In Eq. 2, cont (non-cont) refer to the fact that the residues are (are not) in contact in a given structure or ensemble of structures (two residues are assumed to be in contact if they have any pair of heavy atoms within 5.5 Å) and the symbols G (non-G) indicate that the contact between I and J exists (does not exist) in the native state. As is evident from Eqs. 1 and 2, the validity of a G-type model requires that, for any configuration of the protein,

(3)

that is, that only the native contact interactions contribute significantly to the effective energy for both native and non-native configurations or ensembles. The original form of the G model (Taketomi et al., 1975) assumes that

(4)

where determines whether residues I and J, which make a contact in the native state, are in contact in the structure under consideration, and  is a constant parameter. In the common implementation of G-type models, the assumption is made that the energy for each residue pair in a given structure is proportional to the number N_IJ of native heavy-atom contacts in that structure; i.e.,

(5)

where  is a proportionality constant determined by a fitting procedure (Muñoz and Eaton, 1999; Shimada et al., 2001).

	NATIVE-STATE ANALYSIS

TOP ABSTRACT INTRODUCTION NATIVE-STATE ANALYSIS NON-NATIVE INTERACTIONS VALIDITY OF G MODELS METHODS REFERENCES

Figure 1 shows the distribution of the interaction energy for the interacting residue pairs. It has a peak close to zero, and the average is 1.2 ± 1.8 kcal/mol for the eight proteins. It is evident that the contact energies cover a wide range and that the use of a single coefficient, as in Eq. 4, is a rough approximation.

View larger version (14K): [in this window] [in a new window]   FIGURE 1   Histogram showing the number of residue pairs interacting with a given residue-residue energy. For clarity, results for only four proteins out of the eight studied are shown; their PDB code is 1aps (black), 2ci2 (red), 1hml (green), 1ten (blue). The arrows represent the average interaction energy between pairs of residues, corresponding to the parameter  for each protein.

Figure 2 a shows a scatter plot of the relationship between E_IJ as calculated for the native structure with EEF1 and the number of contacts between residues I and J. Data for four proteins are included in the figure (see caption); the lines represent the best (least-squares) fit for individual proteins and for all the proteins simultaneously. Although a qualitative relationship is evident, there is considerable scatter in the distribution.

View larger version (23K): [in this window] [in a new window]   FIGURE 2   Comparison of EEF1 and G-type energies. (a) Energy between pairs of non-neighboring residues as a function of the number of the heavy atom contacts (with a cutoff of 5.5 Å) between them in the experimental structure. (b) Effective residue energy as a function of the numbers of all (heavy) atom contacts made by that residue. The continuous lines in (a) and (b) are the least-squares fit for the various proteins (thin colored lines) and for all the proteins together (black heavy line). Results are shown for the set of proteins used in Fig. 1.

Table 1 shows the calculated total energies obtained for the native states of eight proteins from the EEF1 energy function in the column headed E(EEF1) (Eq. 1). The next column, E^cont(EFF1), shows the result obtained with the contact approximation, E (Eq. 3), which is equivalent to the exact G-type energy for the native state. E^cont(EFF1) = E^Go(EFF1) is clearly a very good approximation to the true native state energy, E(EFF1). The energies obtained with Eq. 5 are in the next columns. Results obtained by fitting a parameter for each protein (^P) and with an average parameter for all the proteins () are shown; the values of ^P and , are given in Table 2; they correspond to the linear least square fits in Fig. 2a. Deviations from E(EEF1) are as large as 60 kcal/mol with ^P; with significantly larger values occur. We can also introduce a parameter, ^P', which is chosen so that Eq. 5 yields the native state (EEF1) energy for each protein; the values are also listed in Table 2. As can be seen, ^P and ^P' are very similar for each protein (within 0.01 kcal/mol). Nevertheless, because the number of residue-residue contacts is in the range 3000-7000 (as listed in Table 1), such small differences lead to large changes in the total energy. This provides a cautionary note on the use of such formulations.

View this table: [in this window] [in a new window]   TABLE 1   EEF1 total, contact and G energies for the experimental structures of eight proteins (see Methods)

View this table: [in this window] [in a new window]   TABLE 2   Values of P and P used to compute the energies given in Table 1

In applications of G-type models to estimate the effect of single-site mutations on the native-state energy (Otzen et al., 1995; Xu et al., 1998; Cota et al., 2000), and for transition states (Vendruscolo et al., 2001), the quantity of interest is not the interaction energy of residue pairs, but rather the interaction energy per residue, E. This is defined as the sum over all residues in contact with a given residue in the native state (i.e., N = _J N) so that, from Eq. 5,

(6)

This energy is to be compared with the energy per residue obtained from EEF1, which is

(7)

Figure 2 b shows a scatter plot of the results for the proteins included in Fig. 2 a. In analogy to Fig. 2 a, best fit lines for the individual proteins (^P) and for all the proteins simultaneously () are shown. The correspondence is significantly better than that in Fig. 2 a with a narrower distribution about the best fit lines. This result is expected from the central limit theorem, which shows that the standard deviation of the average (or sum) is smaller than the standard deviation of individual variables. The total energies and parameters ^P for the various proteins obtained from the fits to Eq. 6 are given in Tables 1 and 2. The values of ^P are very similar to ^P and to ^P' and, in most cases, are closer than ^P to ^P', the parameters that fit the total energy. This is in accord with the fact that E^Go(^P) is a better approximation to E(EEF1) than is E^Go(^P) in most cases (see Table 1). This is an important result because it provides a justification for the use of Eq. 6 in the analysis of protein stability and transition state structures.

	NON-NATIVE INTERACTIONS

TOP ABSTRACT INTRODUCTION NATIVE-STATE ANALYSIS NON-NATIVE INTERACTIONS VALIDITY OF G MODELS METHODS REFERENCES

The determination of contribution of non-native contacts along folding pathways is important for an evaluation of G-type models. Because the transition states play an essential role in protein folding, we consider them first and then examine other portions of the energy surface.

Transition states

In Table 3, we show the results for the transition state ensembles (TSE) determined by constraining the calculated  values to be equal to the experimental ones, using molecular dynamics and the EEF1 potential (see Methods). The values in Table 3 are obtained by considering a single representative structure of the TSE; corresponding results are obtained if averages over the TSE are made. The structures in the ensembles have a root mean square deviation (RMSD) in the range 4-6 Å from the native state. The first two columns give the EEF1 energy, E(EEF1), used as a reference, and the contact energy calculated with EEF1, E^cont(EEF1). As can be seen, the agreement is as good as it is for the native state (Table 1). The energies calculated using only the native contacts, E^Go(EEF1) are, in all cases, less negative than the true energies. The non-native contribution, E^non-Go(EEF1), also given in the Table, is in the range of 69 to 122 kcal/mol. This shows that non-native contacts contribute significantly to stabilizing the transition state. Table 3 also lists as E^Go(^P) and E^non-Go(^P), the corresponding results obtained using Eq. 5 with ^P, the best fit of the energy parameter to the native state for each protein. There are significant deviations of E^Go(^P) from E^Go(EEF1), in addition to the errors in the latter. The deviations are both positive and negative (between 32 and +60 kcal/mol); in comparison to E^cont(EEF1), the differences are all positive, as for E^Go(EEF1).

Thus, use of the G model, whether in the form of E^Go(EEF1) or E^Go(^P), results in a deeper well for the native state relative to the transition state than do the actual energy values; e.g., for procarboxypeptidase A2 (1aye), the EEF1 transition state energy is 50.5 kcal/mol above the native state, whereas it is calculated to be 124.4 kcal/mol and 92.9 kcal/mol with E^Go(EEF1) and E^Go(^P), respectively. The same is also true relative to the denatured state, in correspondence with the unrealistically deep funnel-like structure of the energy surface obtained with the G-type potential, already mentioned in the Introduction.

Thermally induced non-native conformations

Table 4 shows results corresponding to those in Table 3 for a set of non-native conformations of CI2 obtained by an unfolding simulation at 450 K, followed by quenching simulation at 300 K (see Methods). The contact approximation is valid for these non-native states, as it is valid for native and transition states. However, the various G-type models have errors that increase with the RMSD from the native state; the error in E^Go(EEF1) and E^Go(^P) are similar. For the largest RMSD analyzed (11.7 Å) the stabilization arising from non-native contacts, E^non-Go(EEF1) and E^non-Go(^P) is larger than that from the native (G-type) contacts. This is because there are more non-native (128) than native (76) contacts, as shown in Table 4.

In some simulations (Shimada et al., 2001), it has been assumed that non-native contacts are repulsive, which leads to faster folding than the standard G models, which include only native contacts. This is not in agreement with the EEF1 decomposition, i.e., the non-native constants are overall attractive, as shown in Tables 3 and 4. The G-type potential, used in the folding study of a three-helix bundle protein (Zhou and Karplus, 1999), varied the non-native interactions over a range that included both attractive and repulsive values. A non-native repulsive interaction somewhat weaker than the corresponding native attractive interaction (a ratio of ~0.4) gave a folding time closest to the experimental value.

	VALIDITY OF G MODELS

TOP ABSTRACT INTRODUCTION NATIVE-STATE ANALYSIS NON-NATIVE INTERACTIONS VALIDITY OF G MODELS METHODS REFERENCES

The G model was proposed in 1975 as an ingenious technical device to make possible computer simulations of protein folding (Taketomi et al., 1975). As such, it has been very successful (Takada, 1999). Now the G model is being used not only as a convenient computational tool but also because it is thought by some to capture what they believe to be an essential element (i.e., a smooth deep funnel) of protein energy "landscapes." But does it really do this? Do smooth deep funnels characterize protein energy surfaces?

A comparison between an all-atom molecular-mechanics effective-energy function with solvent shielding and the G-type potentials in current use provides a test of these fundamental questions. It is expected that other effective potentials that are pairwise decomposable would give similar results. The practical impossibility of doing this type of analysis with explicit solvent models should be noted. The results obtained show that the most commonly used G-type model is only an approximate description of the protein potential energy surfaces for the native state, the transition state, and along folding trajectories. Nevertheless, some global features obtained with G-type models are likely to be meaningful. Of particular importance is the demonstration that individual residue energies (rather than residue-residue interaction energies) are rather well described by G-type models. This explains the correlations observed for single-site mutations in the native state (Otzen et al., 1995; Xu et al., 1998; Cota et al., 2000) and provides a justification for the determination of the coarse-grained structure of transition-state ensembles based on such models (Vendruscolo et al., 2001; Paci et al., 2002a). G-type models are most accurate for transition states that are relatively close in structure to the native state, as they appear to be for many fast-folding small proteins (Li and Daggett, 1994; Vendruscolo et al., 2001). For other portions of the potential energy surface, particularly collapsed misfolded states, the non-native contacts neglected in G-type models make important contribution and can lead to significant distortions of the potential energy surface. The present analysis also demonstrates that G-type models result in a much deeper well for the native state than do more realistic potentials, due primarily to the neglect of stabilizing non-native contacts. This is an essential element in the simple, fast-folding behavior obtained in molecular dynamics and Monte Carlo calculations based on G-type models. Thus, just the elements of the G-type models that make them so attractive for folding simulations introduce errors that have to be considered in evaluating the significance of the quantitative folding result for proteins obtained from such model calculations. Moreover, the neglect of attractive non-native interactions makes it impossible to use G-type potentials in the study of misfolding (e.g., the production of fibrils, which appear to be important in certain diseases (Dobson, 1999, 2001)).

	METHODS

TOP ABSTRACT INTRODUCTION NATIVE-STATE ANALYSIS NON-NATIVE INTERACTIONS VALIDITY OF G MODELS METHODS REFERENCES

We use a molecular mechanics potential energy function (EEF1) for the atoms with an implicit solvent term. EEF1 is based on the CHARMM19 polar hydrogen representation (Neria et al., 1996) with a Gaussian model for solvation (Lazaridis and Karplus, 1999). The function, called EEF1, has been used in a variety of applications concerned with the protein folding reaction (Lazaridis and Karplus, 1999), including the high-temperature unfolding of the protein CI2 (Lazaridis and Karplus, 1997), where good agreement was obtained with simulations that used an explicit representation of the solvent (Li and Daggett, 1994).

The effective energy, E^EEF1(R), of a protein with conformation R includes the protein internal energy and the solvation free energy. Both can be written as a sum over all residue pairs. Details are given in Lazaridis and Karplus (1999).

Proteins and conformations used for analysis

Eight proteins were used in the analysis. They are acylphosphatase (PDB entry1aps (Pastore et al., 1992)), chymotrypsin inhibitor 2 or CI2 (PDB entry 2ci2) (McPhalen and James, 1987), -spectrin SRC 3 domain (Blanco et al., 1997) (PDB entry 1aey), the third fibronectin type III repeat from tenascin (Leahy et al., 1992) (PDB entry 1ten), -LA (Ren et al., 1993) (PDB entry 1hml), procarboxypeptidase A2 (Garcia-Saez et al., 1997) (PDB entry 1aye), an immunoglobulin-like modules from titin I-band (Improta et al., 1996) (PDB entry 1tit), the cell-cycle regulatory protein p13suc1, SUC1 (Endicott et al., 1995). The experimental structure was, in all cases, minimized for 200 steepest descent steps to eliminate bad contacts.

Several types of non-native structures of interest for the understanding of protein folding and unfolding were examined. Transition state ensembles were obtained using an approach based on experimental  values to bias the trajectory (Vendruscolo et al., 2001; Paci et al., 2002a) toward conformations where the fraction of native contacts equals the experimental  value for those residues for which such value has been measured. For CI2, high-temperature unfolded states were obtained by increasing the temperature of the Nosé-Hoover thermostat to 450 K during the simulation of over 1 ns or longer. Collapsed configurations were generated from the high temperature conformations by decreasing the temperature to 300 K over 200-ps trajectories.