The inherent structure landscape of a protein

Density of States for the Inherent Structures

Molecular dynamics simulations show that the model has the expected qualitative properties of a good folder: at low temperature it folds to the global minimum of the potential energy, whereas at high temperature it denaturates. If the model is cooled down quickly from a high-temperature unfolded state to a temperature well below its folding temperature (T _f), the folding process is, however, very slow, and may never reach a shape close to the native shape within the time scale accessible in a simulation. This finding suggests that this model protein exhibits the glassy properties observed in actual proteins. The folding transition can be detected from the temperature variation of the specific heat C _v. As shown in Fig. 1, a characteristic peak of C _v appears at a temperature that we define as T _f. The ISL is determined by the basins of the potential energy surface, labeled by an index α, which are defined as the set of conformations that are connected to the same local energy minimum, called the inherent structure, which has an energy e _α.

The rapid change of the shape of the curve showing the probability density of inherent structures P _IS versus their energy e _α in the vicinity of T _f confirms this interpretation: below T _f, P _IS(e _α, T) is higher for the low values of e _α. At T _f it shows two maxima, one for low energies corresponding to folded states, and one for high-energy, unfolded states. Above T_f, P _IS(e _α, T) is dominated by the high-energy, unfolded states. To characterize the ISL we need not only the values of the energies e _α but also their density of states Ω _IS(e _α). A realistic protein has a huge number of inherent states so that Ω _IS(e _α) cannot be determined from a systematic search. Instead it can be deduced from the sampling of the inherent structures that gives P _{IS .} In the basin of attraction of inherent structure α, the potential energy of the protein can be written as V(r) = e _α + ΔV(r), where r designates all of the degrees of freedom of the protein. Let us assume that the inherent structures can be split into a set of discrete states α₀, α₁, …, α _K (which includes at least one state, the ground state) and a continuum of higher-energy states. With this notation, the configurational contribution of the partition function of the protein is where e _max is the highest inherent structure energy. Moreover we have introduced which is the vibrational free energy in the basin of attraction B(α) of the inherent structure α. The probability to be in the basin of attraction of a discrete inherent structure is and, for the continuum of inherent structures, the probability density of inherent structures is If we can assume that the free energy F _v(α, T) does not depend on the inherent structure, in Z the term exp[−βF _v(T)] can be factorized so that p _αi and P _IS simplify into with Z_IS can be viewed as an inherent structure partition function. It can be expressed in terms of the probability p _α0(T) that the protein is in the basin of attraction of the ground state, which is henceforth denoted as p ₀(T). This allows us to rewrite P _IS(e _α, T) as if we chose the ground-state energy as the reference state of the energies (e _α0 = 0). Eq. 7 provides a method to compute Ω_IS(e _α) from the probability density P _IS(e _α, T) deduced from the sampling of molecular dynamics trajectories. The validity of the method relies on the assumption that the vibrational free energy in a basin of attraction F _v(α, T) does not vary significantly from a basin to another. For harmonic basins of attraction, we have where the values ω_q are the frequencies of the vibrational modes of the protein when it is in the configuration corresponding to the inherent structure α. Of course, these frequencies change individually from one inherent state to another, but they keep the same order of magnitude, and the sum has an averaging effect that keeps F _v(α, T) only weakly dependent of the conformation.

This conjecture is verified by our calculations because Ω _IS(e _α) deduced from P _IS(e _α, T) at different temperatures gives the same result in all of the ranges of inherent structure energies that are accessed with a sufficient probability to allow a correct sampling at the temperature considered. Fig. 2shows the density of states Ω _IS(e _α) computed from the sampling of molecular dynamics trajectories at three temperatures. Fig. 2's results are remarkably similar, showing that the assumption that F _v(α, T) can be eliminated in the calculation is a good approximation; however, the results are obviously not identical. Fig. 2 Right, deduced from the calculation at the lowest temperature, shows some values of e _α where the density of states appears to vanish, whereas it does not vanish if Ω _IS(e _α) is computed from a trajectory at higher T. Such inconsistency is particularly noticeable for high values of e _α, and it points out the limit of the numerical calculation. As shown from Eq. 7, the probability to be in the basin of attraction of an inherent structure decreases as exp(−e _α/k _B T). Therefore, if we use a molecular dynamics trajectory at a low temperature to sample the basins of attraction of the inherent structures, it may happen that we completely miss some inherent structures because of insufficient sampling. This outcome is even more likely if the basin of attraction of a given inherent structure is narrow because, besides the temperature factor, its small size also reduces the chance to sample it properly. Studies performed at higher temperatures reduce this problem and T = T _f is particularly appropriate because we expect a molecular dynamics trajectory to sample both structures corresponding to folded and unfolded states. Once Ω _IS(e _α) has been determined from Eq. 7 at a temperature where the sampling is particularly efficient, Z _IS(T) can be computed at any temperature with Eq. 6.

Fig. 3displays Ω _IS(e _α) in a large range of energies. It shows an isolated minimum energy state, followed by a few states very close in energy and then a large gap, of the order of k _B T _f, before the next state. A plot of the structure shows that the ground state and states immediately above it in energy have very similar shapes. Thus the few states near the minimum can be viewed as conformational substates (12) of the ground state. The gap between these states and the next states indicates that the model protein is a good folder (3), as observed in molecular dynamics simulations.

For larger energies Ω _IS(e _α) shows an exponential growth, with two different slopes depending on the energy range:

with

where Ω₀ and A are positive constants, and T _s = 1.1 T _f and T _u = 0.9 T _f. It is interesting to notice that a similar exponential dependence of the density of metastable state was found for spin glasses (13) although the energy dependence of g(e _α) is more complex for such systems.

Thermodynamics in the Inherent Structure Space

The density of states in the ISL can be used to analyze the folding of the protein. Substituting the expression of Ω _IS(e _α) into Eq. 7 gives

But Eq. 11 is also the expression that gives the behavior of P _IS(e _α, T) at low temperature (T < T _u) because, in this temperature range, the equations giving P _IS(e _{α ,} T) show that low-energy inherent structures dominate. On the contrary, Eq. 12 should be used at high temperature (T > T _s), because in this range high-energy inherent structures dominate. As shown in Fig. 1, in the range T _u < T ≈ T _f < T _s none of the energy ranges dominate, and P _IS(e _α, T) has two humps, one of them at low e _α corresponding to folded structures, the other one at high e _α, corresponding to unfolded structures. Therefore, when we take into account the expression of Ω _IS(e _α), we find that the shape of P _IS(e _{α ,} T) changes qualitatively in a narrow temperature range around T _f. This picture is consistent with a first-order transition, T _u < T < T _s being the coexistence region, so that the two temperatures exhibited by the analysis of Ω _IS(e _α) appear as the limit of the spinodals. These characteristics are consistent to the behavior of specific heat C _v. As is indicated in Fig. 1, T _s and T _u are close to the temperatures at the start and the end of the peak for C _v.

As noticed above, the exponential growth for the density of state Ω _IS of the inherent structures is not a special feature of the protein model but was also obtained in spin glasses (13). However, in the present case, the exponent g(e _α) is connected to the folding transition, which is specific to proteins. We also observed a similar behavior in a standard Gō model, without interaction through dihedral angles.

The density of inherent structures energies and the associated partition function Z _IS(T) are not only useful for characterizing the folding transition. They can also tell us a lot about the thermodynamics of the protein. We can build thermodynamics based only on the configurations of the inherent structures by the usual scheme of statistical mechanics,

where U _IS, F _IS, S _IS, or C _IS are the internal energy, Helmholtz free energy, entropy, or specific heat, respectively, deduced only from the thermal fluctuations between conformations, so that they can be called conformational energy, conformational free energy, etc. They are shown in Fig. 4.

The conformational specific heat is remarkably similar to the full specific heat of Fig. 1. The detailed analysis of Fig. 4 Middle, which displays the vibrational specific heat per residue C _W/N _r = (C _v − C _IS)/N _r and conformational specific heat C _IS, shows that the vibrational contribution of the specific heat is very close to the value 3 k _B/2 expected for harmonic oscillators. This finding suggests that the thermalization in the basin of attraction of each inherent structure allows a good separation between a vibrational part and a conformational part of the specific heat. When comparing the values of C _W and C _IS one should notice an essential difference between the two. C _W is a vibrational contribution, which comes from all of the residues and, as it is well approximated by the harmonic value, indicates that the contributions of all vibrational degrees of freedom are independent, so that C _W/N _r has an actual physical meaning. On the contrary, the structural fluctuations entering in C _IS may be either a flip of a residue, for instance, and thus can be very local, or a slow structural change extending to a significant part of the protein, such as the motion of the whole α-helix. This is why the quantity C _IS/N _r would not have a physical meaning because a protein is not a self-similar structure. A meaningful comparison must compare C _IS with C _W/N _r: structural fluctuation will start to play a role in the physics of the protein when C _IS > C _W/N _r. As shown in Fig. 4, C _IS starts to rise significantly ≈0.3 T _f and quickly takes over for T > 0.4 T _f, whereas the fluctuations are dominated by the vibrational contribution in the low temperature range.

The behaviors of the conformational entropy and energy (Fig. 4 Bottom) show a similar rise clearly visible for T > 0.4 T _f. The increase of entropy that appears to accelerate and become exponential for T > 0.4 T _f (see Fig. 4 Bottom Inset) can be understood as coming from the structural variety observed in the thermal fluctuations. This idea would be consistent with the behavior of C _IS and C _W.

Discussion

The results of the thermodynamics of the inherent structure states shed an interesting light on the current debate about the existence and observation of a protein dynamical transition (14). For our model protein two temperatures appear to play a particular role. The lowest one is 0.05 T _f, above which energy and entropy show a rise (Fig. 4). With the resolution of our calculation,^§

The results from thermodynamics in the inherent structure space are not affected by a possible lack of ergodicity at low temperature because they are established from a density of states deduced from simulations around the T _f. Their accuracy is determined by the quality of the sampling of the phase space in the vicinity of the T _f.

<0.05 T _f, fluctuations among inherent states seem completely suppressed, so that this temperature could be viewed as the Kauzmann temperature for this model (15). One should, however, notice that Fig. 4 shows a gradual decrease of the entropy ≈0.05 T _f, suggesting that, as previously noticed for another off-lattice protein model (16), a true thermodynamics glass transition does not exist for a finite system such as a protein. The other characteristic temperature is 0.4 T _f. It should correspond to a qualitative change in some observations on the protein because, around this temperature, conformational changes among the different inherent structures start to dominate the fluctuations, as shown by the specific heat. Thus a glassy behavior could be expected when the temperature is lowered to values in the range 0.05 T _f < T < 0.4 T _f or below.

However, as far as we deal with equilibrium aspects, the characteristics of a glass phase are not detected in our study because they only show up in out-of-equilibrium situations. From the viewpoint of equilibrium thermodynamics, we cannot find a characteristic change caused by the “glass transition.” In the temperature range 0.05 T _f < T < 0.4 T _f, in which we could expect the glass transition, the specific heat C _v only shows a smooth evolution without any drop when temperature is lowered. In the thermodynamic sense there is no transition. However, the picture could become different if nonequilibrium aspects were probed. For instance, when cooling the protein rapidly from high temperature, a drop of C _v might be observed and a transition temperature identified in the range 0.05 T _f < T < 0.4 T _f, but it should depend on the cooling rate.

Our results raise a question about the nature of the fluctuations that we observe: do they correspond to the α or β relaxations observed in proteins (14)? Experiments (17, 18) and molecular dynamics simulations in which the temperature of the solvent and the protein are constrained to be different (19) show that the solvent mobility is an important factor in determining the large-scale fluctuations of a protein. The α relaxation is “slaved” to the solvent fluctuations, which is not the case for the glassy aspects of our results because the solvent only enters into our study through its contribution to the interaction energies not through its dynamics because the model uses effective potentials to include hydrophobic effects. The fluctuations must be attributed to the protein energy landscape, which suggests that the fluctuations that are observed in our simulations correspond to the β relaxation, i.e., to fluctuations that do not involve large conformational changes that could be hindered by a frozen solvent in an actual experiment. This makes sense if we think of the peculiarities of the Gō model, which strongly favors the native structure by the selection of the interactions that are introduced in the model. It is likely that such a bias toward the native shape prevents such a simplified model from exhibiting an α relaxation.

In the studies of glassy systems with inherent structures, it is customary to introduce the configurational entropy

In the context of a thermodynamic description based on the inherent structures, s _c can be viewed as the microcanonical entropy, whereas S _IS(T) defined by Eq. 15 is the canonical entropy. Therefore, in the thermodynamic limit, we could expect that s _c(ẽ_α) = S _IS(T), where the most probable value ẽ _α of the inherent structure energy would also tend toward its mean value 〈e _α〉. However, this is not true for a protein, because of the finiteness and heterogeneity of the system, but also because, below T _f, the most probable value of e _α is the ground state of the protein. In this regime the density of states Ω _IS(e _α) is dominated by a few isolated states, and one cannot build thermodynamics on the configurational entropy, whereas S _IS(T) defined by Eq. 15 has the expected properties of a thermodynamic entropy, such as dS _IS(T)/dU _IS(T) = 1/T.

In this article we have shown that the ISL of a protein can be analyzed because its density of states can be deduced from the probability density of inherent states extracted from molecular dynamics simulations. The ISL contains considerably less information than the free-energy landscape of the protein, and this is what makes it accessible. At the expense of much longer calculations, the ISL appears to be accessible for an all-atom model of a simple protein because our scheme to obtain it does not require the full determination of every local minimum because it takes advantage of the usual ability of statistical physics to provide very accurate results from a simplified description of the system. An important aspect is that, in this framework, all thermodynamics properties are deduced from the density of inherent states that is obtained from molecular dynamics simulations in the vicinity of T _f to sample folded and unfolded states. As a result the difficult issue of the ergodicity of simulations at very low temperature can be avoided.

When a complicated model, such as an all-atom description is considered, the landscape becomes immense and one may be concerned by the validity of the sampling of the phase space to get the inherent structure density of states. But it should be noted that, even for a frustrated Gō model, the landscape is rugged and complex. For our protein model, it is multifunneled. Kinetic studies of folding detect the existence of a second folding funnel, leading to a kinetic trap (unpublished work). The density of inherent states in the kinetic trap can be obtained and it is very low. Therefore, although the second funnel plays a significant role in the kinetics of the folding, it does not show up in the thermodynamics of the protein. This is an interesting aspect of the canonical sampling that we use. It introduces some self-consistency in the method because states that are important for the reduced thermodynamics are also those that are more likely to be sampled. Another test of the validity of this approach is provided by the possibility to get Ω _IS(e _α) by sampling at different temperatures in the vicinity of T _f. They should lead to the same structure for the density of states. This is why, despite the obvious difficulty of carrying long molecular dynamics trajectories at different temperatures for a complex protein model, the thermodynamics in the ISL seems reachable in this case.

The Gō model of a real protein is sufficient to allow the investigation of the properties of the ISL. It shows that, despite the loss of information with respect to the full free energy landscape, the ISL still contains a lot of data on the protein. It can be used to study the folding transition because P _IS exhibits a qualitative change across T _f and can even detect the associated spinodals. The density of inherent states shows an exponential scaling, which can be related to the folding transition. This finding suggests that it could be a general property of proteins, and it will be interesting to test this hypothesis on other protein models, including all-atom descriptions. The ISL can also be used to analyze the thermodynamics of the protein. It cannot observe the vibrational contributions, but detects the most important features of proteins, their conformational fluctuations that are responsible for their function.

Materials and Methods

The potential energy of our Gōmodel of protein G is

where N _r is the number of residue, and b _i, θ _i, and φ _i are the ith bond length, ith bond angle, and ith dihedral angle, respectively. r _ij is the distance between ith and jth residues. The index 0 denotes the parameters determined from the native structure. Such a model includes local and nonlocal interactions. The local interactions are the bonds that define the C^α chain through stiff elastic potentials describing the covalent bonds, angular potentials between adjacent bonds, and dihedral potentials. The natural bond length or bond angles are determined from the native structure. The nonlocal bonds connect C^α carbons that are not adjacent along the peptide chain, but nevertheless close in space in the folded geometry of the protein. Such carbons are said to form a “native contact” and interact with a Lennard-Jones potential in the model. They are determined from the experimental structure of the protein in its folded state: two C^α carbons are said to form a native contact if the distance between them or between atoms of the side chains attached to them is less than a critical distance, here chosen as 5.5 Å. A repulsive interaction is introduced between all C^α carbons that do not form a native contact to describe the steric repulsion. Such a Gō model is specifically designed to ensure that the system reaches the correct folded structure and these models are generally good folders, i.e., their folding transition is sharp and leads to a well defined structure, because they are minimally frustrated. A comparison of the folding of an off-lattice Gō model and a frustrated model based on hydrophobic (B), hydrophilic (P), and neutral (N) residues (BPN model) showed that the Gō model leads to a sharp transition where collapse and folding occur simultaneously, whereas the highly frustrated model shows a very broad collapse and even at the lowest temperature the number of contacts is far from maximal (20), pointing out that a realistic model of protein should be minimally frustrated (2). In our model a slight frustration is introduced (10) by assuming that the dihedral angles have the ideal values 45^° or 135^° of the α-helices or β-sheets instead of the angles deduced from the native structure. The native structure of protein G has been obtained from the Protein Data Bank (21). The dimensionless parameters of the model are K _b = 200.0, K _a = 40.0, K _d = 0.3, ε = 0.18, and C = 4.0. The experimental structure lies in the basin of the global minimum of the model but does not correspond exactly to the ground-state structure of our model because of the frustration introduced on the dihedral angles.

Thermalized phase space trajectories are generated for this model, using underdamped Langevin simulations, where the masses of all of the residues are assumed to be equal. Our goal here is not to explore out-of-equilibrium states, but rather to analyze the equilibrium properties of a protein. Therefore, to avoid long nonequilibrium transients that occur on cooling, except for the general tests of the model properties, all of our simulations started from the ground state and heated up to the temperature of interest. To collect statistics over a given equilibrium trajectory, a typical simulation included 3 × 10⁸ steps after the desired temperature was reached. For each temperature 20 simulations were performed with different sets of random numbers in the Langevin equations.

Equilibrium trajectories were used, on one hand, to analyze the thermodynamic properties of the model, and, on the other hand, to characterize its ISL. The basins of the potential energy surface were determined by quenching 6 × 10⁴ instantaneous configurations separated by 10⁵ steps, i.e., calculating a steepest-descent path from these points of the trajectory. A similar algorithm had been used earlier for the BPN model of a four-helix bundle protein (22) to detect states that could be long-lived intermediates in the folding.