Energy landscape and dynamics of an HP lattice model of proteins —The role of anisotropy

Understanding proteins is a daunting challenge because of the inherent complexity of the problem. Yet, there has been major progress in the development of a theoretical framework by borrowing ideas from spin glass theory and transition state theory [1–4]. Simplified lattice models [5–9] have proved to be invaluable in obtaining exact results that have proved to be a guide for interpreting experiments. Recently, a Go model [10,11] was studied for a linear chain of twelve beads on a square lattice [12]. The Go model is generally meant to be a simplified way of incorporating minimal frustration and maximal compatibility. The key results of the exact numerical studies of the energy landscape and the folding dynamics of the twelve-bead model were at odds with conventional beliefs. Specifically, they were surprising on several counts: a) the energy landscape was characterized by kinetic partitioning [13,14], b) the dynamics of folding was not the standard two-state behavior of folded and unfolded conformations and instead was controlled by trap states, c) the φ values did not represent the native like nature of the transition states [15,16] but rather depended sensitively on the conformations of the native state, the dominant trap state, as well as the numerous transition states between the trap state and the native state. Of course, the Go model is merely a means of encoding the native state with artificial interactions between the native state contacts. The model studied in ref. [12] also lends itself to a conventional HP model with two kinds of beads hydrophobic (H) and polar (P). Our aim here is to study such an HP model exactly and discuss what kinds of modification might make the model more protein-like.

We consider a 2-dimensional variant of the HP model [9] in which the beads are either hydrophobic (H) or polar (P). In order to reduce degeneracies, the H-H, H-P, and P-P contacts are assigned energies of −1, −0.2, and −0.1 respectively. There are 2¹² possible sequences and some of these have a non-degenerate ground state conformation that coincides with the previously studied native state state NAT of the Go-like model and which is shown in the inset in the top panel of fig. 1. Among the most stable among these sequences is PHHPHHHPPHHH corresponding to the native state energy $E=-5.1$ and the folding transition temperature (at which the equilibrium native state occupancy is $\frac{1}{2}$) is $T_f=0.091$. The density of states of this system is shown in the top panel of fig. 1.

The dynamics of folding is carried out numerically following the exact procedure presented ref. [12]. This procedure involves a numerical integration of the Master equation [17,18] for 15037 conformations. The system exhibits 403 local energy minima and 117 extended minima (the latter comprise 2407 conformations) resulting in a very rough energy landscape compared to the Go-like system. The extended minima are sets of equal-energy states which are mutually accessible to each other but getting to any other conformation requires a supply of energy. This system does not fold at T_f over a time scale of $100000\ \tau$ when the native state is made a sink of the dynamics. Here, τ denotes the characteristic microscopic time in the Master equation. In contrast, for the Go-like system of ref. [12] folding was achieved within $2000\ \tau$. The lack of folding is not surprising because many of the local energy minima have exit barriers which are substantially larger than T_f. Many of the minima have a barrier of 0.2, including the lowest local minimum, which is at −4.4. Two minima with $E=-4.2$ have a barrier of 2 and 1. Such minima become very long-lived traps.

A basic summary of the results is that the HP model has an energy landscape that is not protein-like in its attributes. Because of low-lying states close to the native state, the temperature at which there is significant native state occupancy is too low to overcome the barriers from the trap states thereby leading to glassy behavior and thwarting rapid folding. We now turn to a discussion of what kinds of generic modifications one can make to the HP model to remedy this problem. What one needs to accomplish is to change the energy landscape to make the native state deeper and the folding transition temperature higher. The dynamics of folding could then improve especially if the folding transition temperature is high enough for any barriers associated with significant trap states to be overcome.

In order to make progress, we look to Nature to provide a clue. The key secondary structures build on a simple attribute of any linear chain molecule. Two chain segments placed parallel to each other are able to sustain interactions along their length whereas two chain segments that are not parallel are limited in the regime in which they are spatially near each other. In a helical arrangement, the successive turns of a helix lie alongside each other. Likewise, adjacent strands of a sheet structure are also alongside each other. This shows that the inherent tube-like anisotropy of any chain molecule prefers parallel placement of neighboring tube segments [19].

This element of the tube picture of proteins is readily incorporated in the tube HP model (THP) [20]. This is accomplished by relating the strength of the contacts to their context by enhancing the attraction between parallel segments of the chain. Specifically, we multiply the contact energies by a factor, which is equal to either 3, 2 or 1, depending on whether the two interacting beads belong to segments which are fully, partially, or not parallel, respectively. (A somewhat similar Go-like scheme based on the degree of nativeness is discussed in ref. [21].) The ground state of the system corresponds to the same conformation NAT found in the Go-like and HP systems but its energy is now −14.2. The resulting energy spectrum is shown in the bottom panel of fig. 1. It comprises 417 local energy minima and 115 extended minima (housing 2360 conformations). These number characteristics are similar to those of the HP model and yet the THP system does fold at its T_f of 0.993 as shown in the inset in the lower panel of fig. 1. The reason for the different kinetic behavior is that the exit barriers are generally small compared to T_f. The largest among these is equal to 1.4 and is associated with the lowest excited local energy minimum at −11.4 —this state is denoted as TRAP in fig. 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The TRAP state controls the folding process. Among 25 conformations of energy −8 there are 10 which are the lowest transition states between NAT and TRAP in the sense that there are T = 0 trajectories, which lead from these states to both NAT and TRAP (not at equal rates).

Even though the THP system results in folding, the time scale to folding is 45% longer than in the Go-like model because of the larger roughness in the energy landscape in the former. When one quenches from the extended states, essentially no probability ends in NAT. The probability gets trapped in many non-native states. The two states that catch the probability the most are shown at the top of fig. 2. Thus the kinetic partitioning into the native state is absent unlike in the Go model. Another way to characterize the landscape connectivity is to begin with every state (i.e. also those which are not extended), one at a time, and monitor whether there is a T = 0 path that can lead to NAT. We find that only 18.7% of all conformations are connected kinetically to NAT in this way, compared to 94% in the Go-like model.

We observe that the THP model improves folding over the HP model but it still comes with a rough energy landscape. An inspection of the location of the energy levels shows that the roughness the landscape stems directly from the adoption of non-zero values of energies of the H-P and P-P contacts this is a convenient choice to ensure that degeneracies of the energy levels are decreased within the standard HP model. However, for the THP model, this step is not necessary. The conformation denoted as NAT constitutes a unique ground state within THP even if the contacts involving P are assigned zero energy. We denote a modified THP model with non-zero contact energies ascribed only to H-H pairs as the THP' model.

The top panel of fig. 3 shows that the energy landscape in THP' is much smoother than for the same sequence considered within the THP model (the bottom panel of fig. 1). The landscape comes with 65 local energy minima instead of 417. The T_f is nearly the same as in the THP model, but the folding process at T_f is faster —the TRAP state, shown in fig. 4, has an exit energy barrier of 1 instead of 1.4. Furthermore, at T = 0, 13.4% of probability flow, starting with open conformations with no contacts at all, lands in NAT which restores the kinetic partitioning found in the Go model [12].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Remarkably, one can smoothen the energy landscape further by considering a symmetric sequence such as the one shown in the bottom panel of fig. 3. It comprises 6 H-beads and it has the same unique ground state which coincides with NAT of the previous cases. but now it corresponds to the energy of ${-}14$. The energy landscape now has just 28 local energy minima and 4 extended minima. T_f is somewhat lower −0.81. The resulting foliding time is about 30% shorter than for the THP' model of the assymmetric sequence. There is now substantial kinetic partitioning as NAT collects 47% of the probability at T = 0 from the open conformations. The trap states, shown in the bottom of fig. 4 are of energy −9, they jointly collect 52% of the probability, and their exit barriers are equal to 1.

In summary, we have carried out numerically exact studies of the energy landscape and the folding dynamics of a simple HP model and the corresponding tube-like HP model for a 12-bead chain on a square lattice. These studies along with results of previous work on a Go-like model with the same ground state lead to surprising results. The standard HP model does not exhibit protein-like characteristics and is unable to fold in measurable time scales. In contrast, the incorporation of the chain-like anisotropy, albeit in an approximate way, leads to more protein-like behavior. We have studied different variants of the tube HP model leading to different degrees of kinetic partitioning to the native state. Strikingly, the THP model with a symmetric sequence yields a more funnel like energy landscape than the Go model. Such a funnel necessarily involves kinetic partitioning involving significant barrier-less routes to the native state. As a corollary, the rate limiting steps in folding necessarily involve escape from trap states and thus the φ values reflect the routes from the trap states to the native state through transition states. Thus all simplified models of proteins studied exactly, be it the Go model or the THP model or its variants, do not exhibit standard two-state behavior and the canonical interpretation of the φ values does not hold. Our results suggest that underneath the seeming simplicity of behavior of small fast folding proteins might lie some unexpected complexity.

The computer resources were financed by the European Regional Development Fund under the Operational Programme Innovative Economy NanoFun POIG.02.02.00-00-025/09. This research has been supported by the Polish National Science Centre Grant No. 2011/01/B/ST3/02190.