Landscape approaches for determining the ensemble of folding transition states: success and failure hinge on the degree of frustration

Abstract

We present a method for determining structural properties of the ensemble of folding transition states from protein simulations. This method relies on thermodynamic quantities (free energies as a function of global reaction coordinates, such as the percentage of native contacts) and not on "kinetic" measurements (rates, transmission coefficients, complete trajectories); consequently, it requires fewer computational resources compared with other approaches, making it more suited to large and complex models. We explain the theoretical framework that underlies this method and use it to clarify the connection between the experimentally determined Phi value, a quantity determined by the ratio of rate and stability changes due to point mutations, and the average structure of the transition state ensemble. To determine the accuracy of this thermodynamic approach, we apply it to minimalist protein models and compare these results with the ones obtained by using the standard experimental procedure for determining Phi values. We show that the accuracy of both methods depends sensitively on the amount of frustration. In particular, the results are similar when applied to models with minimal amounts of frustration, characteristic of rapid-folding, single-domain globular proteins.

Figure 1

Approximate schematic of a Φ value, which holds for proteins that are weakly frustrated. The solid curve is a schematic free energy profile for a wild-type protein; the dashed curve is for a suitable mutant. The free energy profile is drawn against a single order parameter for folding. For many small proteins, many simple global structural measures of nativeness may be used for this free energy projection. ΔΔG⁰ and ΔΔG^‡ are, respectively, the change in native state stability and change in activation free energy upon mutation. A Φ value near 1 suggests that the local environment of the mutated residue is native-like in the transition state; a Φ value near 0 suggests that the local environment of the mutated residue is unfolded-like in the transition state. This interpretation of Φ values becomes less valid as the frustration increases.

Figure 2

The free energy profile G(Q) shown at T_f for the three models studied in this paper and one of their mutants. The three models have the same structure but different potentials: a Gō-like potential where every bead is a different type, and only beads adjacent in the native structure are attractive; a 3LC sequence ABABBBCBACBABABACACBACAACAB, where contacts of identical type have energy of −3 units and contacts of different type have energy of −1 units (arbitrary energy units are used to have a folding temperature and glass temperature of order unity); and a 2LC sequence ABABBBBBABBABABAAABBAAAAAAB. The mutant (M11) is from the 3LC system, in which the native interaction between beads 6 and 13 (numbered starting from the lower left corner of the structure shown in the inset) is reduced in energy from −3 to −1. All mutants show two-state behavior. (Inset) The native 3 × 3 × 3 structure for all of the different wild-type and mutant sequences studied in this paper.

Figure 3

A comparison of the kinetic and free energy perturbation methods for inferring folding transition state structure. Each panel shows the comparison for a different sequence (and Hamiltonian): A is the Gō-like sequence, B is the 3LC sequence, and C is the 2LC sequence. The Gō-like sequence has the least amount of energetic frustration, and C the most. The native structure and potential for these sequences are shown in Fig. 2. Mutants are made by weakening specific nonbonded interactions between beads that are adjacent in the native structure. The Φ values for these contacts are then computed by the standard experimental procedure (ordinate) and by a free energy perturbation technique (abscissa). Agreement is good (normalized correlations 0.86 and 0.84) for the models in A and B, which have energetic frustration less than or comparable to small, fast-folding globular proteins. Because the Gō sequence in A has no energetic frustration, the heterogeneity of the Φ values is mostly determined by topological factors attributable to a combination of the polymeric nature of the chain and the structure of the native state. More frustrated sequences, such as in C, show no agreement (normalized correlation −0.49) between the two methods and out-of-range Φ values, which suggests that the assumption of a Kramer's type of rate with a fixed rate prefactor is not valid. In the experimental method (ordinate), the folding rate of the wild type (k_wt) and the mutant (k_mut) as well as the change in native stability under mutation are measured and used to compute a Φ value as Φ −RT ln(k_mut/k_wt)/ΔΔG⁰. These Φ values should be similar to the measure of ΔΔG⁰/ΔΔG^‡, the ratio of the change in the folding activation free energy to the change in native stability, when the assumption that folding follows a Kramer's type of equation with a fixed rate prefactor is valid. In the free energy perturbation method, we determine the free energy as a function of a folding reaction coordinate—in this instance, Q, the fraction of formed native nonbonded interactions. The barrier height is defined as the difference in free energy between the highest free energy point along Q between 5/28 and 23/28 and the free energy of states, with Q < 16/28. The Φ value is then computed directly from ΔΔG⁰/ΔΔG^‡ by taking Q < 16/28 as the unfolded conformations. For both methods, the unfolded state is defined as all conformations with Q < 16/28. Error bars show 68% confidence limits calculated from 1,000 bootstrapping simulations.

Figure 4

A comparison of the 28 bond Φ values—produced by mutating a single native interaction, decreasing its energy from −3 to −1—from the model with a Gō-like potential and a 3LC potential. (The sequence and native structure are shown in Fig. 2.) The agreement shows that, for sequences with reduced frustration, the native structure (topology) has a large role in determining the Φ values. Details of the potential interactions may not be as critical. Φ values are determined here from the rates by using Eq. 1.

Figure 5

(Upper) A comparison of the Φ values computed for the Gō-like model from the free energy perturbation formula (based on Eq. 4) by using two different formulas for estimating the activation free energy ΔΔG^‡. The abscissa shows the Φ values computed by assuming a fixed transition state location at Q = 16/28; the ordinate shows the Φ values computed by assuming that the activation free energy is equal to the variation between the unfolded free energy minimum and the maximum barrier point. The solid line is a least squares fit constrained to pass through the origin with a slope of 1.07. (Lower) The same plot but shown for the 3LC model. The solid line has a slope of 1.38. These two plots demonstrate that the less frustrated models have less sensitivity of their Φ values to the assumed position of the barrier. The overestimation of the Φ values computed from free energy perturbation in the 3LC model is apparent from the lower panel. Because the top of the free energy barrier is very broad, variations in the reconfigurational diffusion coefficient with Q and the existence of fundamental motions that allow jumps of several units in Q can shift the actual barrier location. In this instance, the actual location is close to Q = 16/28, which is at the lower Q position of the barrier.

Figure 6

A comparison of the Φ values computed by measuring changes in the folding rates and by averaging over the thermally weighted population of conformations with transmission coefficients between 0.4 and 0.6. Comparison is done for the minimally frustrated Gō sequence and structure shown in Fig. 2. The two results have a normalized correlation coefficient of 0.73—agreement is thus of comparable quality to using our method with a single order parameter Q but at a much greater computational cost. Ten-thousand states with Q values in the range 6/28 to 27/28 were sampled in equilibrium at ≈100,000 step intervals. Three-hundred and ninety states from this sample had transmission coefficients between 0.4 and 0.6, defined as the fraction of folding simulations that, when started in a given conformation, find the native conformation before any conformation with Q ≤ 7 (as determined from 100 independent folding simulations). This set of structures was then used as a putative transition state ensemble for computing Φ values via the free energy perturbation equation.