2.1 Evaluating NMA's predictive power for real proteins

In some cases, the protein architecture – steric interactions arising from spatial distribution of density, coordination number and chain connectivity – is so important that the structure itself may be predictive of other possible functional conformations and even of pathways between them. This occurs more frequently with quaternary motions for large multidomain proteins or protein oligomers. Hemoglobin and kinases like AdK may serve as examples (Arora & Brooks, Reference Arora and Brooks2007; Karplus & McCammon, Reference Karplus and McCammon1983; Whitford et al. Reference Whitford, Miyashita, Levy and Onuchic2007). Hence, when a large protein can be fruitfully thought of as comprising several rigid bodies that are loosely connected (Karplus & McCammon, Reference Karplus and McCammon1983; Tama & Brooks, Reference Tama and Brooks2006), the large-scale slow motions are likely dominated by relative displacements of these ‘rigid bodies’ (subunits or domains) as a whole with respect to each other (Fig. 2). Atoms inside the rigid bodies are tightly packed, characterized by high coordination numbers. These cases are identified by a low number of contacts in the interface between the subunits, which are changing in such motions, compared with the nearly unchanged contact maps within the rigid bodies. Consequently, the energy of contacts is not changing significantly in the transition, and therefore the effect of conformational entropy change is expected to be of significant importance.

Fig. 2.

Some special types of protein molecular architecture allow a mechanistic view on the protein's motions. In this schematic representation, the protein can be thought of as a machine that has two rigid parts connected by a hinge, with motion around that hinge being a single degree of freedom.

To clarify the last statement, let us imagine a rod with one end attached to the origin that can rotate around the z-axis at a fixed angle, θ (Fig. 3). If a slow change of this angle is associated with a conformational motion, then the rotational entropy will be affected too, since it is proportional to sin θ. Analogously, such ‘rigid-bodies-on-hinges’ motions in various protein complexes may be partially driven by changes in the protein conformational entropy.

Fig. 3.

Entropy contained in a degree of freedom depends on the spatial Jacobean. In this figure, the rotational entropy of the rod becomes larger with larger opening angle, because the size of space available to the rod's end is proportional to the sine of this angle.

One successful example of the ‘rigid-bodies-on-hinges’ technique application was a study of hemoglobin, where the first normal mode turns out to be relevant to going from the tense (T) to the relaxed (R2), oxygen-bound state (Xu et al. Reference Xu, Tobi and Bahar2003). If one starts from the T state and applies a small perturbation along the slowest mode of motion, as defined by the anisotropic network model (ANM), then the root mean square distance (RMSD) to the R2 state will decrease as the magnitude of this step along the first mode grows, from 3·5 Å (RMSD between T and R2) down to 2·5 Å. In other words, the motion along the first mode from the T state approaches the R2 state. A step in the other direction along the mode takes the protein farther from R2, and there is no way of ascertaining the motion direction, without knowing the R2 structure in advance. But even leaving a choice of two structures is indicative of the method's predictive power. Analyzing the weights with which the coordinates of different residues enter the first mode, authors found that the least mobile atoms are on the interface between α₁ and β₂ subunits, hinting that they mediate a hinge-like motion. The authors suggested that the transition from T to R2 is favored entropically. As the contacts between the subunits (hydrogen bonds and salt-bridges), responsible for maintaining the T state, weaken upon oxygen binding, the hemoglobin transitions to the R2 state primarily because of increased conformational entropy. In the analogy above, this motion is similar to the increase of the angle θ between the rod and the z-axis, along with the corresponding increase in entropy. This view is corroborated by the fact that the α₂β₂ unit is rotated by more than 15° upon the conformational change from T to R2. Interestingly, when the authors repeated the procedure backwards, starting from R2, then the first mode did not lead to the T state, in either direction. This result is not surprising, as the backwards transition cannot also be driven by an increase in conformational entropy. On the other hand, because of the lack of a detailed mechanism, the picture of entropically driven transition still has to be regarded critically.

In the above example, elastic network modes not only serve as a good basis set for the motion description but are also able to predict the other functional conformation, also suggesting a specific pathway for the transition. However, it is not clear how often other proteins would produce such favorable outcomes, even if they consist of nearly rigid domains connected by hinges. It is likely that the role played by the architecture is extreme in this example, based on the specific topology of steric interactions and chain connectivity. Furthermore, the hemoglobin transition is well described by the first mode only, which is quite unusual. In other proteins, even if the architecture leads to a transition driven by conformational entropy, a combination of several low-frequency modes typically contributes to the allosteric transition pathway. This diminishes the predictive power of the technique with respect to finding the alternative conformation. It is relatively simple to trace a single mode and produce a series of structures, but it is less practical to cover a space of several dimensions where each point is a possible structure.

Nonetheless, the success of describing allostery with NMA is not unique for hemoglobin. Another similar example is the T→R transition in the chaperonin protein GroEL (Tehver et al. Reference Tehver, Chen and Thirumalai2009). A quantitative measure for the success of a particular mode in the description of the transition (Alexandrov et al. Reference Alexandrov, Lehnert, Echols, Milburn, Engelman and Gerstein2005; Marques & Sanejouand, Reference Marques and Sanejouand1995; Tehver et al. Reference Tehver, Chen and Thirumalai2009) is the overlap defined as

(1)

where x _i^T and x _i^R are the atomic coordinates in states R and T, respectively, and a _M is the Mth mode. This quantity equals the absolute value of the cosine of the angle between a _M and x ^R−x ^T, covering a range from 0 to 1. When the cosine is equal to 1, at zero angle, it is indicative of complete coincidence of the mode and the difference between T and R structures. For the T→R transition in a single subunit of GroEL, the overlap for the lowest frequency mode was found to be 0·83 (Tehver et al. Reference Tehver, Chen and Thirumalai2009). Thus, the first mode could have been used to predict the R structure if the latter had not been known.

In addition to detailed studies of allosteric motion for particular proteins, the robustness of the harmonic motion approach has been systematically tested on larger protein datasets (Alexandrov et al. Reference Alexandrov, Lehnert, Echols, Milburn, Engelman and Gerstein2005; Tama & Sanejouand, Reference Tama and Sanejouand2001). An analysis by Tama et al. used a set of 20 proteins for which both ‘open’ and ‘closed’ forms are known (Tama & Sanejouand, Reference Tama and Sanejouand2001), calculating overlaps of the normal modes using Eq. (1). The mode that has the maximal overlap is considered ‘the most involved’ in the conformational change. In 7 out of 20 cases, this was the lowest frequency mode. Thus, in the majority of cases, the use of the first normal mode for the prediction of the other allosteric conformation would not have been successful. Overlap in those seven proteins varies from 0·3 to 0·72 with an average of 0·51. The index of ‘the most involved mode’ was not very small for all the 20 proteins. For seryl-tRNA synthetase, for instance, the most involved mode had index 18. The overlap of the most involved mode with conformational difference (Eq. 1) was ranging from 0·3 to 0·83, with an average of 0·56. In a different study, Alexandrov et al. (Reference Alexandrov, Lehnert, Echols, Milburn, Engelman and Gerstein2005) filtered all possible pairs of ~33 000 sequences from Protein Data Bank (PDB) for sufficient sequence identity (>95%), structural dissimilarity (>1·5 Å RMSD) and then, through CATH (Class, Architecture, Topology, Homologous superfamily) database, identified 377 pairs that should have a conformational transition between them. In this study, overlap was defined in a slightly different way: Eq. (1) (calculating a cosine of angle between atom displacement along the normal Mth mode and the difference in position in two real structures) was applied for each atom individually. In addition, as an integral measure, the root mean square of these individual overlaps was computed. When these calculations were performed without using additional criteria (such as only looking at atoms that move the most), the average angle between the vectors was only marginally different from that between random vectors. However, this study likely included many smaller, single-domain proteins that possess only tertiary motions, which are ill-defined by the protein architecture.

NMA provides a better mode transition overlap when the initial reference conformation is the open form (Tama & Sanejouand, Reference Tama and Sanejouand2001). This is not surprising, as the shape and lack of connections are more pronounced in those forms. Additionally, when a conformational transition is more collective, overlaps also tend to be higher (Tama & Sanejouand, Reference Tama and Sanejouand2001). All these facts indicate that harmonic mode-based techniques are most successful when applied to the problems similar to the hemoglobin or GroEL examples (Tehver et al. Reference Tehver, Chen and Thirumalai2009; Xu et al. Reference Xu, Tobi and Bahar2003) – to study quaternary motions of large multidomain proteins or oligomers, preferably with a single hinge region. However, even in these cases, truly successful examples have been rare.

A large number of NMA studies use the reproducibility of the B-factors – the experimentally determined values indicating the amplitudes of atomic fluctuations in a crystallized protein – as a way to validate the technique, or for parameterizing the model (Bahar & Rader, Reference Bahar and Rader2005; Eyal et al. Reference Eyal, Chennubhotla, Yang and Bahar2007; Hayward & Go, Reference Hayward and Go1995; Kondrashov et al. Reference Kondrashov, Cui and Phillips2006; Riccardi et al. Reference Riccardi, Cui and Phillips2009; Soheilifard et al. Reference Soheilifard, Makarov and Rodin2008; Tehver et al. Reference Tehver, Chen and Thirumalai2009; Zheng, Reference Zheng2008). However, individual atomic fluctuations (characterized by B-factors) are mostly defined by the atom's local environment, such as atomic coordination numbers, and in densely packed parts of the protein these local environments will not change significantly regardless of the motions of these parts. Therefore, reproducing B-factors is not an impressive test of the validity of applying elastic network models to a particular protein. However, the fluctuations obtained by NMA are sensitive to the cut-off radius of the network model, so B-factors can be profitably used for choosing the spring constants or the cut-off radius. Caution is still needed: in protein crystals, the local atomic environment is different than in solutions, because of interactions with neighboring molecules in the lattice; therefore, crystallographic B-factors may noticeably differ from the atomic fluctuations in the solution (Kundu et al. Reference Kundu, Melton, Sorensen and Phillips2002; Soheilifard et al. Reference Soheilifard, Makarov and Rodin2008), while nuclear magnetic resonance (NMR) results show better correspondence with fluctuations in elastic network models (Yang et al. Reference Yang, Eyal, Chennubhotla, Jee, Gronenborn and Bahar2007). B-factors are determined by contributions from motions on all scales, both spatially localized and collective, while the main physical idea of NMA is to study slow large-scale motions, which are largely defined by protein architecture. In other words, the springs are representative of the architecture, not the actual energy landscape above the ångström scale, as we discuss further in the section.

2.2 Adiabatic normal modes, double wells and cracking

We now consider a hypothetical protein, having two native conformations A and B with a barrier separating them. Motivated by hemoglobin and other allosteric proteins and enzymes, we are interested in the dynamics of the transition between these states. Problems of this kind are frequently addressed by NMA techniques. Since each of the two states is a free energy minimum, it is natural to depict an approximate one-dimensional free energy surface comprising two quadratic wells with a barrier between them. In reality, the landscape is multidimensional and allows multiple transition paths, but to simplify the problem one may start with two conjugated parabolae, in an approach similar to that used in the electron transfer theory (Marcus & Sutin, Reference Marcus and Sutin1985, Miyashita et al. Reference Miyashita, Onuchic and Wolynes2003).

One must always keep in mind that the main reason why the normal modes may serve as suitable reaction coordinates for large-scale protein motions is the fact that these motions are largely determined by the constraints imposed by the protein architecture. While it is physically reasonable to assume that the native conformation (typically determined by X-ray crystallography) is a minimum of the energy of contacts, and expand these energies to the second order, this harmonic approximation will likely hold in a very narrow vicinity of either structure A or B. Therefore, the NMA springs connecting the atoms are mainly reflecting the architecture of the conformation rather than the actual energy landscape.

If architecture is indeed paramount, NMA may still guide the motion further from the initial state (A) to the final state (B) even when the harmonic approximation fails. It is possible then to adiabatically modulate the normal modes while constructing the transition path, that is, slightly change the structure, as normal modes suggest, and then calculate new normal modes, repeating this process until the final state is reached. This approach has been used with adenylate kinase, where the putative transition between open and closed forms was iteratively followed by NMA steps (Miyashita et al. Reference Miyashita, Onuchic and Wolynes2003). First, normal modes in state A and their overlaps with the transition displacement vector were calculated. The overlap is a scalar product of a normal mode 3N-dimensional vector and the difference between target and initial structures (3N-dimensional vector whose components are the differences between the atoms' coordinates in states B and A). Then a small step along the highest overlap mode (or a combination of a few highest overlap modes) was taken. Subsequently, the normal modes were recalculated for the new state (conformation), and the new step forward was taken. This procedure generates a sequence of conformations that may reflect the actual transition path, and a reaction coordinate along that path, λ, related to the number of steps along adiabatically changing normal modes. Naturally, the procedure can be also repeated in the opposite direction, starting from state B and going along the modes that overlap with the A–B structural difference (Miyashita et al. Reference Miyashita, Onuchic and Wolynes2003).

The procedure described above does not guarantee that similar forward and backward trajectories will be generated, nor is it clear how to estimate their correspondence to the actual transition path. Nonetheless, all the conformations from the iterative sequence of steps can be mapped on a traditional reaction coordinate, such as RMSD (root mean square deviation of atomic coordinates) or Q (Plotkin et al. Reference Plotkin, Wang and Wolynes1997) (the fraction of shared contacts with the initial structure), or two coordinates, e.g. Q _A and Q _B (Okazaki et al. Reference Okazaki, Koga, Takada, Onuchic and Wolynes2006; Wu et al. Reference Wu, Zhuravlev and Papoian2008), or RMSD(A) and RMSD(B), like Miyashita et al. did. To eliminate the second coordinate, authors established a correspondence between RMSD(B) and RMSD(A) for each point of the transition path based on the path length (Miyashita et al. Reference Miyashita, Onuchic and Wolynes2003). It would be interesting to build a 2D free energy surface in terms Q _A and Q _B from all-atom MD (Wu et al. Reference Wu, Zhuravlev and Papoian2008), and compute a dominant path on this surface to compare it with the forward and backward paths obtained with the iterative NMA approach.

Another subtle point with this scheme is that in addition to obtaining the transition path itself, one often would like to know the free energy profile along that path to predict the kinetic rates. To address this problem, the authors build on an approach analogous to electron transfer theory and assume two conjugated parabolic potentials. If higher accuracy is required, a free energy profile may also be computed by performing an umbrella sampling along λ using a detailed MD force field. This should be feasible by reversing the correspondence between λ and RMSD. In addition to mapping the adiabatic normal mode pathways in two dimensions, as mentioned above, this approach might lead to a logical extension of the iterative NMA technique, providing higher resolution details.

The next question is how far can normal mode techniques be pushed while still providing meaningful results. Real energy landscapes are almost always poorly known, and parabolic wells are their simplest approximation. Energy landscapes of globular proteins are funnel-like (Bryngelson et al. Reference Bryngelson, Onuchic, Socci and Wolynes1995; Dill & Chan, Reference Dill and Chan1997; Dobson et al. Reference Dobson, Sali and Karplus1998; Frauenfelder et al. Reference Frauenfelder, Alberding, Ansari, Braunstein, Cowen, Hong, Iben, Johnson and Luck1990; Onuchic et al. Reference Onuchic, Wolynes, Luthey-Schulten and Socci1995; Papoian et al. Reference Papoian, DeGrado and Klein2003a; Plotkin & Onuchic, Reference Plotkin and Onuchic2002; Yue & Dill, Reference Yue and Dill1995), that is, defined by the native contacts. As these contacts are stretched, the free energy grows, and it is plausible that this growth is quadratic with respect to the atomic coordinate displacements. Furthermore, the spring constants can be tuned to correspond to the actual funnel slope. To save this simple representation beyond the point where the harmonic approximation breaks down, Miyashita et al. suggested that the main reason for its violation is that the parabola is too steep, and it becomes more favorable for the chain to locally unfold, or crack the contacts, than to continue stretching them (Miyashita et al. Reference Miyashita, Onuchic and Wolynes2003, Reference Miyashita, Wolynes and Onuchic2005). To incorporate this possibility into the model, the authors calculated the strain energy in each conformation along the adiabatic normal modes path as the strain of the springs of the initial structure (A or B). The strain energy turned out to be localized in sequence, thus supporting the cracking hypothesis. The authors assigned an entropy gain to each ‘unfolded’ residue and compared that to the energy in the case where the residue remained ‘folded’ (this energy consists of contact energy and strain energy accumulated over adiabatic deformation). Local unfolding changes the free energy surface from quadratic to locally linear, and the same change occurs for the dependence of the barrier height on the free energy difference between the two states: without the cracking it is quadratic, similar to Marcus theory for the electron transfer, while it becomes linear when the cracking is allowed (Miyashita et al. Reference Miyashita, Wolynes and Onuchic2005) (Fig. 4). This fact may be used to quantitatively test the cracking model.

Fig. 4.

Cracking is the breaking of some of the native contacts and partial unfolding. If there are two allosteric states, then in harmonic approximation, there are two quadratic surfaces, and the rise up of their slopes may be less favorable than partial unfolding. In this case, the molecule might crack and transition to the other surface at a lower free energy.

The idea of conjugating potentials for two different protein states is quite general and may be applied to other models as well. In the case of Gō-like force fields, the conjugation is also related to local unfolding (Okazaki et al. Reference Okazaki, Koga, Takada, Onuchic and Wolynes2006). We elaborate on this issue further in this paper.

2.3 Phonons

On very short (picosecond) timescales, proteins do behave like rigid bodies and exhibit phonon dynamics based on atomic vibrations. This range is experimentally accessible through infrared (IR) spectroscopy, and a number of computational techniques permit interpretation of the experimental results. The problem of the extreme locality of the harmonic approximation remains at these timescales. Instantaneous NMA counteracts this problem by generating a number of solvent cages and performing NMA for each one of them (La Nave et al. Reference La Nave, Scala, Starr, Sciortino and Stanley2000; Schmitz & Tavan, Reference Schmitz and Tavan2004; Schrader et al. Reference Schrader, Schreier, Cordes, Koller, Babitzki, Denschlag, Renner, Lweneck, Dong, Moroder, Tavan and Zinth2007; Schulz et al. Reference Schulz, Krishnan, Daidone and Smith2009). Energy flow in proteins may be described in a harmonic approximation through quantum heat transfer operators (Leitner, Reference Leitner2008, Reference Leitner2009), or by just observing classical heat flow through a frozen protein globule in MD simulation (Ota & Agard, Reference Ota and Agard2005). These techniques shed light on communication pathways within a protein, which we review in the next chapter.

3.1 Pre-existing equilibrium

Most of the current literature dedicated to describing the mechanisms of allostery is based on the ‘pre-existing equilibrium’ approach, related to the MWC model, as discussed earlier. In the last two decades, it has become widely accepted that the native state consists of an ensemble of conformations on an energy landscape rather than a single static state, and the idea of the pre-existing equilibrium among these states is very natural. Likewise, from this perspective, it is more appropriate to talk about the change of the native ensemble, rather than a conformational change (Benkovic et al. Reference Benkovic, Hammes and Hammes-Schiffer2008).

Let us consider the joint energy landscape of a protein and its ligand (Fig. 5, upper panel). With no ligand bound, the protein vigorously explores the left part of this landscape, spending most of the time in conformation A (which is experimentally observed when no ligand is bound), only occasionally visiting conformation B (which is experimentally observed when the ligand is bound), in the state denoted as β on the joint landscape. During one of these visits to the β state, a ligand binds, locking that conformation, so the system moves to state B on the joint landscape. Importantly, the conformational step is much slower than the binding step. Pre-existence of an equilibrium assumes multiple transitions occurring between A and β on microsecond to millisecond timescales. Diffusion of the ligand to the binding site, once the ligand is within an ångström from the binding interface, is on the scale of nanoseconds. For this reason, the MWC mechanism is often called a population shift mechanism, since it can be envisioned as sudden deepening of the minimum at β (becoming B) with simultaneous shallowing of the minimum at A as it goes to α (Arora & Brooks, Reference Arora and Brooks2007; Okazaki & Takada, Reference Okazaki and Takada2008). On the joint landscape, we can think of population shift as the system jumping from the unligated cross-section of the landscape (with A and β) to the ligated cross-section (with B and α). Note that there does not have to be a pronounced minimum at α, since it is not biologically significant.

Fig. 5.

Upper left: a surface with four minima is the simplest schematic representation of the joint ligand-and-protein energy landscape. Upper middle: in the MWC scheme, the protein is in equilibrium between states A and β in the absence of the ligand and upon ligand binding goes to state B, which is the same conformation of the protein as β but with a ligand. The first step is much slower (μs–ms) than the second (ns). More likely, the MWC transition would follow a less extreme pathway, similar to the one denoted as ω. A minimum α does not have to exist, but it facilitates interpretation of the MWC picture as a ‘population shift’. In the KNF scenario, binding starts first, and there is no timescale separation between binding and conformational change. This may happen along a pathway similar to the one denoted as κ. Upper right: the previous statements can be summarized in a phase diagram: the lower right zone of the plot hosts KNF pathways (the start of binding precedes conformational change) and the upper left zone hosts the MWC pathways (conformational change precedes binding). Lower left: to avoid considering the complicated joint energy landscape, the MWC picture may be visualized as a diabatic transfer between two energy surfaces: one with a ligand and one with no ligand, since the switch between the surfaces happens much faster than visitation of the state β on the unligated surface. Such an interpretation also suggests that the transfer will possibly occur earlier than the visit of state β is completed, as it is too high in free energy. Lower right: in the KNF scheme, 1D interpretation is a free energy profile along the transition pathway. If this pathway on a joint landscape is enveloped by a valley, the transition can be thought of as a gradual adiabatic relocation of the minimum in which the protein continues to reside. Note: many ingredients that are used to construct this figure have been employed long ago. However, the authors of this manuscript have introduced some new elements, and the overall synthesis has many new aspects. It has to be seen whether the theoretical description of allostery proposed in this review will be validated by detailed simulations and experimental work.

Figure 5 (lower left) shows another interpretation of an MWC-type allosteric transition that does not consider the joint energy landscape of a protein and its ligand, but only that of the protein. It is analogous to diabatic electron transfer between two surfaces. The surfaces in the case of the MWC switch are the energy landscapes of the protein molecule in unligated and ligated states (they are largely the same as the cross-sections of the joint landscape discussed above). The unligated surface has another minimum at β, but the same conformation has a much lower energy when the ligand is bound, on the other surface. Hence, if the ligand becomes available while the protein happens to explore the neighborhood of β, the protein almost instantaneously jumps to the other surface. Such an immediate jump is called a diabatic transfer. Depicting the MWC mechanism in this way suggests that the diabatic transfer may happen prior to the protein chain fully reaching the β state, since two diabatic surfaces are likely to cross at lower energies before that point. Therefore, if the ligand has approached, the dominant transition pathway may be something like the one denoted as ω or bold arrow (Fig. 5) rather than actually two-step A→β, β→B.

Indeed, it seems to be the limiting case of the MWC picture, when conformation B is visited with extreme precision in the absence of the ligand, having correct conformation of the backbone and local side-chain orientations at the interface primed for binding the ligand. Experimental evidence to date does not possess sufficient structural resolution to ascertain if this is a likely scenario. We consider a number of cases below in detail.

The concept of a pre-existing equilibrium started to develop vigorously after an indication from MD simulations that myoglobin shows the largest conformational fluctuations in the regions that are also the most different in the X-ray structures for allosteric states (Frauenfelder et al. Reference Frauenfelder, Parak and Young1988). Early experimental evidence has been reviewed by Kern & Zuiderweg (Reference Kern and Zuiderweg2003). It includes the observation of exchange between active and inactive states in unphosphorylated NtrC by NMR spectroscopy (Volkman et al. Reference Volkman, Lipson, Wemmer and Kern2001) (allosteric transition occurs on phosphorylation), a transient loop disengagement in the R-state of fructose-1,6-bisphosphatase (the disengaged loop corresponds to the T-state) revealed by X-ray and fluorescence (Nelson et al. Reference Nelson, Iancu, Choe, Honzatko and Fromm2000), the observation that the conformation of CO-ligated hemoglobin (defined by relative orientation of the α₁/β₂ subunits) shown by NMR is exactly halfway between T and R crystal structures interpreted as dynamic average between them (Lukin et al. Reference Lukin, Kontaxis, Simplaceanu, Yuan, Bax and Ho2003).

More direct evidence appeared with time. The small angle X-ray scattering curve for wild-type aspartate transcarbamylase changes abruptly at about 60°C, indicating the transition, while the same curve for a mutant changes gradually and non-monotonically, suggesting that it represents a dynamic average of two states, whose relative populations shift slightly with temperature variation (Fetler et al. Reference Fetler, Kantrowitz and Vachette2007). The kinetic curves for glucose binding by human glucokinase are best explained by the kinetic scheme that suggests pre-existing equilibrium of states with higher and lower affinity for glucose (Antoine et al. Reference Antoine, Boutin and Ferry2009; Kim et al. Reference Kim, Kalinowski and Marcinkeviciene2007). Similar schemes also explain fluorescence anisotropy curves for cyclic AMP (cAMP) receptor protein of Escherichia coli (Youn et al. Reference Youn, Koh and Roberts2008).

An MD study of acetylcholinesterase classified 89 PDB structures from five different species and grouped them according to χ₁ and χ₂ angles of Trp279, which is the key residue of the peripheral anionic site that may interact allosterically with the catalytic anionic site and may show conformational change upon binding (Xu et al. Reference Xu, Colletier, Jiang, Silman, Sussman and Weik2008). The authors find five islands on the (χ₁, χ₂) plot. Some of the groups of PDB structures are seen on the islands, that is, corresponding rotamers are observed in the simulation, so these rotamers are in pre-existing equilibrium. However, the conformational state of the rest of the molecule was not included in the clustering analysis.

Evidence for a pre-existing equilibrium has been also observed in vivo (Ayers et al. Reference Ayers, Nedrow, Gillilan and Noy2007). FABP4 is a protein that performs transport from the cytosol to the cell nucleus. To enter the nucleus it must have a nuclear localization signal (NLS) conspicuously exposed. In FABP4, NLS consists of three residues, and the conformation where NLS is at work, allowing the entrance to the nucleus, is locked by a regulating ligand. However, even without the regulating ligand, some accumulation of FABP4 was observed in the nucleus, suggesting that the corresponding allosteric conformation is visited even in the unligated state (Ayers et al. Reference Ayers, Nedrow, Gillilan and Noy2007).

It has been shown that a mutation in cyclophilin A can invert relative populations of the pre-existing substates (Fraser et al. Reference Fraser, Clarkson, Degnan, Erion, Kern and Alber2009). This mutation also slows conformational dynamics and decreases catalytic turnover by two orders of magnitude. Further spectroscopic investigation suggests that enzymatic activity decrease is linked to slowed conformational dynamics. As in the last example, NMR spectroscopy has revealed global conformational rearrangements in G-proteins that are thought to be the reason allowing the protein to interact with diverse ligands (Smrcka et al. Reference Smrcka, Kichik, Tarrago, Burroughs, Park, Itoga, Stern, Willardson and Giralt2010).

In summary, a small number of works have directly supported the ‘pre-existing equilibrium’ hypothesis. In many cases, however, only local conformational rearrangements were considered, where it is not possible to positively claim that the other allosteric conformation is visited in its entirety. Nonetheless, the picture of equilibrated ensemble of states is very natural in the energy landscape paradigm. Hence, some studies use this language of ensemble equilibrium and population shift without any pointers as to which allosteric mechanism is relevant to their particular case (Martinez et al. Reference Martinez, Gohon, Corringer, Tribet, Merola, Changeux and Popot2002; Malmendal et al. Reference Malmendal, Evenas, Forsen and Akke1999). This language has even led to a suggestion that allostery is a purely thermodynamic phenomenon (Tsai et al. Reference Tsai, Sol and Nussinov2009). As previously mentioned, it is true that there is always an ensemble rather than a static conformation, but there is also a biological timescale, which may not accommodate the visitation of the other allosteric state in the absence of the ligand. Therefore, even if these cases are fundamentally the same from the statistical mechanics viewpoint, where finite barrier heights are irrelevant, they may be very different from the biological, or practical, perspective. In other words, the pre-existing equilibrium picture may not hold for all proteins: some would have landscapes with higher barriers, which cannot be crossed on the timescales corresponding to the observed kinetic rates, and therefore binding is needed to physically induce the transition. Visualizing this situation within the energy landscape paradigm, one would observe continuous relocation of the minimum from conformation A to conformation B during the ligand binding process rather than pulling down of a pre-existing local minimum at B. On the other hand, ‘relocation of the minimum’ is a rather inaccurate term. Because of timescales involved in KNF, the ‘relocating minimum’ may be very local. On the joint ligand-and-protein landscape, the dominant KNF pathway would probably be similar to the one denoted as κ (Fig. 5). It should be occurring fast, on a nanosecond to microsecond timescale: by definition, the key step in this mechanism is binding (nanosecond timescale), which may be slowed down to some extent by the need for conformational rearrangement. Therefore, if one defines the allosteric transition coordinate along the path κ, the equilibrium transverse to it will only be on the timescales faster than 10⁻⁸ s. Nonetheless, since there is no timescale separation between binding and conformational change, it is convenient to think of the transition as adiabatic, with a free energy profile along the dominant transition path similar to the one shown in Fig. 5 (lower right).

Thus, the joint ligand-and-protein energy landscape can be roughly partitioned into two areas, corresponding to the two mechanisms, reminiscent of a phase diagram (Fig. 5, upper right). There is no strict boundary between the zones, but largely if the ligation starts earlier than conformational change (the systems move to the right before moving upwards), that corresponds to the KNF mechanism and the opposite corresponds to the MWC scenario.

Structurally, a range of intermediate situations between pure MWC (A→β→B) and KNF may also be possible or even likely. During the dynamics without a ligand, conformations other than A may be visited and they may resemble B without being exactly the same (Nelson et al. Reference Nelson, Iancu, Choe, Honzatko and Fromm2000; Xu et al. Reference Xu, Colletier, Jiang, Silman, Sussman and Weik2008). Resemblance to the B conformation could be complete, including orientations of the side chains in the whole molecule, or only at the level of similarity of the backbones to just a local reorientation at the active site, so that an enzyme could accommodate a ligand – and, subsequently, the induced fit via reaction front would take over from there. On the joint landscape in Fig. 5 (upper right), these scenarios would correspond to the pathways starting in the MWC zone (upward) and, then turning closer to the KNF zone, maybe even crossing into it. Again, from this perspective the ω pathway seems to be more natural than the extreme A→β→B case.

Returning to the picture of the gradually changing protein energy landscape upon binding in the KNF scheme, the path κ in Fig. 5 corresponds to the case when there is a well-defined reaction coordinate for the transition. In that case, there is a dominant pathway (something like κ), and a free profile along that pathway should be similar to the one given in Fig. 5 (lower right). The landscape has to be intricately built to provide a dominant pathway. Possibly, in some cases, it is partially formed by frustration localization near the binding interface (Ferreiro et al. Reference Ferreiro, Hegler, Komives and Wolynes2007). A single reaction coordinate is not the only option. The total timescale range in the KNF scenario (from the transition time (10⁻⁷ s–10⁻⁶ s) to the atomic vibrations (10⁻¹⁵ s)) is smaller than the one in the MWC scenario, so there is less space for a timescale gap that could help us to single out slower, important coordinates. It is possible, therefore, that one may need a considerable fraction of the conformational coordinates to properly describe multiple transition pathways.

Finally, intrinsically disordered proteins, which fold upon binding of a ligand, present a special case. Their folding landscape is often weakly funneled (Papoian, Reference Papoian2008), but in the state without a ligand, the backbone conformational entropy wins over the folded active structure, masking it among a myriad of other conformations, despite the fact that this structure can be quite low in energy.

3.2 Communication pathways

The main evidence supporting the reaction front picture is the experimentally found existence of ‘communication pathways’ across the bodies of allosteric proteins. For some allosteric protein families, which are evolutionarily related, variations of residues in certain positions turned out to be correlated with variations of residues in other positions, sometimes, quite distant in space (in 3D structure). One of the most studied examples is the PDZ domain family. In 1999, Lockless and Ranganathan approached the problem of interaction of distant residues from the evolutionary perspective (Lockless & Ranganathan, Reference Lockless and Ranganathan1999). Double-mutant cycle to measure real thermodynamic interactions between the residues had existed for some time, but evolutionary conservation of such interactions in different homologues of the same enzyme would suggest that the interaction is required for function. Based on a large set of sequences a statistical ‘free energy’ was introduced. Mapping of this statistical mutation free energy (based on 274 PDZ domains) on a 3D structure of PDZ domain resulted in a continuous path consisting of residues with high statistical mutation free energy and going across the whole globule (Lockless & Ranganathan, Reference Lockless and Ranganathan1999). The authors also introduced a statistical free energy of interaction between two positions, based on statistics from the same set of 274 sequences. These free energies turned out to be correlated with thermodynamic free energies of residue interactions obtained by mutant cycle analysis for one residue position of one of the PDZ domains. This coincidence points to the assumption that the identified pathway is indeed based on energetic couplings, although without revealing the coupling mechanism, which would be important for enzyme functioning. Subsequent double mutant cycle works on a PDZ domain, however, have challenged the claim of correlation between statistical and thermodynamic free energies (Chi et al. Reference Chi, Elfstrom, Shi, Snall, Engstrom and Jemth2008).

In 2004, the actual dynamics of the residues along the communication pathway in the PDZ domain were elucidated by NMR (Fuentes et al. Reference Fuentes, Der and Lee2004). In the second PDZ domain from human tyrosine phosphatase 1E (hPTP1E), certain residues (connected by van der Waals interactions, forming thus a contiguous ‘path’) change their ps–ns side-chain dynamics upon binding of the ligand. This change propagates up to ~11 Å away from the ligand. The observed ‘path’ correlated with the evolutionarily conserved path reported by Lockless & Ranganathan (Reference Lockless and Ranganathan1999).

The PDZ domain is not the only case where a network of evolutionarily conserved residues is found (Benkovic & Hammes-Schiffer, Reference Benkovic and Hammes-Schiffer2003; Suel et al. Reference Suel, Lockless, Wall and Ranganathan2003; Zheng et al. Reference Zheng, Brooks, Doniach and Thirumalai2005). Mutations in residues distant in space from the binding site can change the enzymatic rate quite significantly, by up to a factor of 200 (Benkovic & Hammes-Schiffer, Reference Benkovic and Hammes-Schiffer2003), suggesting some form of communication between these distant residues and the enzyme active site. Whether these communication pathways have been optimized by evolution is an interesting question.

In many works, the communication is envisioned in a mechanistic fashion, using phrases like ‘binding energy is not randomly dissipated, but is directed’ or ‘mechanical deformation’ that propagates along the pathway. It is true that sterically connected atoms allow propagation of mechanical or heat waves (Leitner, Reference Leitner2008, Reference Leitner2009; Ota & Agard, Reference Ota and Agard2005), but we recall that, first, the timescale for complete energy dissipation is on the order of fs–ps and, second, binding energy is actually a free energy. Thus, the latter, being directed along a pathway, may, for example, be stored, in side-chain entropy of the residues forming the pathway. The timescale of communication propagation is not known at this time. Mechanistic or quantum heat wave pictures are one possibility, but this mechanism is unlikely if the communication occurs on times longer than several picoseconds. We hypothesize another possible mechanism of front propagation: that the residues along the pathway are frustrated (Bryngelson et al. Reference Bryngelson, Onuchic, Socci and Wolynes1995; Ferreiro et al. Reference Ferreiro, Hegler, Komives and Wolynes2007; Onuchic & Wolynes, Reference Onuchic and Wolynes2004; Onuchic et al. Reference Onuchic, Wolynes, Luthey-Schulten and Socci1995; Whitford et al. Reference Whitford, Gosavi and Onuchic2008a) (and thus predisposed to a change) and the binding of a ligand induces a ‘flipping’ of the pathway resides. An experimentally observed allosteric switching mechanism in caspase-1 may be reinterpreted using this idea (Datta et al. Reference Datta, Scheer, Romanowski and Wells2008).

Caspase-1 is a dimeric enzyme with two cooperatively working active sites on opposite sides of the dimer and an allosteric site at the center. In each subunit, there is a contiguous network of 21 hydrogen bonds formed by nine side chains and a salt bridge in the on state (when the active site is occupied). In the off state, there is a different network, with only 12 hydrogen bonds (only one bond is present in both networks). Thus, the hydrogen-bonding network rearranges upon binding, while maintaining a site-to-site connection (Fig. 6). These residues are well-conserved evolutionarily, and their mutations by alanine scanning show that four specific residues are important for enzymatic activity. Mutation of any of the two residues forming the salt bridge reduces the catalytic activity by 100–200-fold. Two other residues showed minor 4–7-fold activity reduction effect upon mutation (Datta et al. Reference Datta, Scheer, Romanowski and Wells2008).

Fig. 6.

The dimer caspase-1 has two binding sites connected by a hydrogen bond network that rearrange upon binding. In the figure, the sites are occupied by ligands (shown in purple). The residues participating in the hydrogen bond network are highlighted in red and blue. Red residues form a continuous path, hinting at the possibility of the communication pathway. Of utmost importance are the two residues on the interface between the subunits that form salt bridges. Their mutation to alanine drops the enzymatic rate by a factor of 200 (Datta et al. Reference Datta, Scheer, Romanowski and Wells2008).

For the PDZ domain in particular, the results of a very interesting computational work may be viewed as an indirect argument for the fast propagation of the allosteric signal. PSD-95, a member of the PDZ domain protein family, was cooled down to 10 K in an MD simulation without solvent, and then the key active site residue (His76) with ligand and one water molecule bound was heated to 300 K (Ota & Agard, Reference Ota and Agard2005). The absence of explicit solvent is required to weakly hold the other atoms in their places and to introduce a distance-dependent dielectric constant. The authors observed the heat transfer through protein and found it to be highly anisotropic and in agreement with pathways obtained by Lockless & Ranganathan (Reference Lockless and Ranganathan1999). The timescale of this thermal diffusion was on the order of a few picoseconds. The fact that heat travels this fast means that there is a densely packed ‘tube’ along the direction of propagation. However, it does not necessarily follow that an allosteric signal propagates as fast, despite following the same pathway. In addition, despite general correspondence, there were differences between pathways obtained in the simulation (Ota & Agard, Reference Ota and Agard2005) and the statistical analysis of sequences (Lockless & Ranganathan, Reference Lockless and Ranganathan1999) (although, according to the authors, it may be attributed to the fact that the simulation was done for a particular sequence, PSD-95, while the statistical analysis was done using 274). A similar approach is possible, where some localities in a molecule are ‘pumped’ by a harmonic perturbation (10–100 ps scale), and other parts are being ‘probed’ for response (Sharp & Skinner, Reference Sharp and Skinner2006).

The idea of a 1D communication pathway that is used for allosteric signal propagation is intriguing, but it requires more direct evidence to become an established paradigm. As we argue next, the results described so far largely correlate with the existence of a high-density region(s) within the protein core. The residues in the evolutionarily conserved contiguous path reside in the high-density core and thus have a high burial free energy, which is picked up by double mutant cycle analysis. Therefore, it is also not surprising that they are more conserved, as the change in the core is energetically expensive. The NMR studies do point to a contiguous pathway of residues dynamically responding to the binding, but it does not directly follow that this is a pathway for allosteric conformational change.

Although the pathway concept appeals to our everyday experiences, a front propagating through 3D medium is most naturally expected to be a 2D front. In reality, of course, it can be a manifold of any dimension lower than three (including fractal) down to zero. A zero-dimensional reaction front corresponds to the communication pathway picture and seems quite exotic from this point of view. In fact, many studies, even before 1999, were starting from seeking correlation in structural or dynamical change in distant residues upon binding (e.g. any allostery, not necessarily a pathway across the protein) (Benkovic & Hammes-Schiffer, Reference Benkovic and Hammes-Schiffer2003; Gianni et al. Reference Gianni, Walma, Arcovito, Calosci, Bellelli, Engstrom, Travaglini-Allocatelli, Brunori, Jemth and Vuister2006; Hilser et al. Reference Hilser, Dowdy, Oas and Freire1998; Lenaerts et al. Reference Lenaerts, Ferkinghoff-Borg, Stricher, Serrano, Schymkowitz and Rousseau2008; Niu et al. Reference Niu, Chen, Zhang, Shen, Shi and Wu2007; Zhuravleva et al. Reference Zhuravleva, Korzhnev, Nolde, Kay, Arseniev, Billeter and Orekhov2007). They would often find that many distant residues ‘respond’ to a residue near the active site, and vice versa. Inside the nascent communication pathway paradigm initiated by Lockless & Ranganathan (Reference Lockless and Ranganathan1999), such results were formulated as the presence of ‘communication networks’ inside a protein molecule. For instance, 14 conservative (i.e. deleting methyl groups in the side chain without altering the stereochemistry of the residue) mutations to the second PDZ domain from the tyrosine phosphatase PTP-BL (PDZ2) led to changes in the rate constants, despite the fact that the deleted side chains were not participating in the interactions with the ligand (Gianni et al. Reference Gianni, Walma, Arcovito, Calosci, Bellelli, Engstrom, Travaglini-Allocatelli, Brunori, Jemth and Vuister2006). Authors state that most of the affected residues form a network of hydrophobic interactions in the core. However, four similar mutations in the third PDZ domain of PSD-95 did not affect the rate constants, so the effect is PDZ2 specific (Gianni et al. Reference Gianni, Walma, Arcovito, Calosci, Bellelli, Engstrom, Travaglini-Allocatelli, Brunori, Jemth and Vuister2006). Residues showing a change in millisecond dynamics, measured by NMR, in AF-6 PDZ domain upon ligand binding were also named a dynamical network (Niu et al. Reference Niu, Chen, Zhang, Shen, Shi and Wu2007). In barnase, 11 residues that are distant from the binding site show a change in 10–100 ps dynamics upon binding of barstar (Zhuravleva et al. Reference Zhuravleva, Korzhnev, Nolde, Kay, Arseniev, Billeter and Orekhov2007).

Thus the paradigm leads to the temptation of representing the complicated distribution of dynamical or structural perturbations inside a protein molecule as a large number of communication pathways. From a MD simulation in the CHARMM force field with some restrictions derived from NMR experiments, Dhulesia et al. report the presence of two networks in the second PDZ domain of human tyrosine phosphatase 1E that consist of several hundreds of pathways (Dhulesia et al. Reference Dhulesia, Gsponer and Vendruscolo2008). With the total number of residues in the PDZ domain being just below 100, it is not clear whether this picture of many intertwined pathways does anything to help envision the reaction front propagation. Assuming that a protein should be envisioned as a spaghetti coil of ‘communication cables’, perturbing one of them would stir the rest, inducing cross-talk across the whole molecule; but the reported pathways consist only of 4–5 residues and most of the residues are involved in many pathways. In this case, it seems more natural to think of the single higher dimensional fractal manifold formed by the perturbations.

Therefore, caution is needed when exploring communication pathways by either experimental or computational means. However tempting, one should not assume their prevalence in proteins before further direct proof is obtained. For instance, a fast and inexpensive computational method for finding the communication pathways was developed and applied to study the power stroke of myosin (Tang et al. Reference Tang, Liao, Dunn, Altman, Spudich and Schmidt2007). The method is based on seeking the shortest contiguous path of physically interacting residues, which are also evolutionarily conserved. This method is guaranteed to find pathways, but further tests will be needed to examine their relevance to allostery. On the other hand, the hydrogen bond network in caspase-1 is indeed a site-to-site pathway, and that pathway is both conserved evolutionarily and important for activity (Datta et al. Reference Datta, Scheer, Romanowski and Wells2008), but there is no evidence that any reaction front travels along it. In fact, the authors of this alanine scanning study discussed their results in the framework of the pre-existing equilibrium mechanism (Datta et al. Reference Datta, Scheer, Romanowski and Wells2008).

The NMR methods provide the most direct evidence of pathways in a PDZ domain; however, they do not show pathways for all proteins, for example, calmodulin (Igumenova et al. Reference Igumenova, Lee and Wand2005) or eglin c (Clarkson et al. Reference Clarkson, Gilmore, Edgell and Lee2006), where some of the mutations invoke contiguous response, while response to others is non-contiguous. Uncoupling between the dynamics of certain residues upon transition (that had been dynamically coupled before transition) may also play an important role (Clarkson et al. Reference Clarkson, Gilmore, Edgell and Lee2006; Lenaerts et al. Reference Lenaerts, Ferkinghoff-Borg, Stricher, Serrano, Schymkowitz and Rousseau2008). We suggest that signal propagation also may be engendered by specific distribution of local frustration across the molecule. A recently developed method allows the measurement of frustration localization (Ferreiro et al. Reference Ferreiro, Hegler, Komives and Wolynes2007) and may help to shed light on this issue. Finally, the ps timescale in the case of mechanical signal propagation over a high-density pathway is too removed from the allosteric timescales, which are usually in the ns–ms range.

We conclude this section by suggesting a picture that synthesizes some of the discussed mechanisms: communications may occur via a front of highly frustrated residues, forming a ‘frustration tube’, with fractal dimensionality likely higher than one, where the tube residue conformations ‘switch’ upon binding of a ligand. This possibly might be the case for the caspase-1 dimer, discussed above (Datta et al. Reference Datta, Scheer, Romanowski and Wells2008). In general, it may be expected that different mechanisms operate to couple distant residues in various allosteric proteins and enzymes. Even if the original communication pathway paradigm does not prove to be ubiquitous, it has spurred very fruitful discussion among the community, bringing under focus the issue of complicated communication networks in proteins. In particular, these questions are raised directly in the real space, as opposed to the picture based on collective coordinates, which had been more common.

4.1 Revealing the landscape topographies

One straightforward, yet computationally demanding approach to constructing a representation of the protein energy landscape is based on isolating conformational states from the protein dynamics (usually, realized as an MD or Monte Carlo simulation), assigning statistical weights to these states, and computing rates for inter-state transitions. This approach avoids the need for selection of a few collective coordinates, where the latter tends to provide a rigid frame for the whole study. For this reason, the choice of collective coordinates is extremely important and, as discussed earlier, can be highly non-trivial. On the other hand, mapping the conformational space onto a network only requires an order parameter that is good locally, providing a criterion for clustering the multitude of recorded structures into network nodes (conformational states). Since these requirements are looser, it is usually easier to find such a criterion. It may be based on analyzing RMSDs (Krivov & Karplus, Reference Krivov and Karplus2004), contact maps derived from positions of C_α atoms (Hori et al. Reference Hori, Chikenji, Berry and Takada2009) or mutual overlaps, Q. The edges of the network may be derived by analyzing the actual transitions in a simulation (Caflisch, Reference Caflisch2006; Khalili & Wales, Reference Khalili and Wales2008; Krivov & Karplus, Reference Krivov and Karplus2002, Reference Krivov and Karplus2004; Krivov et al. Reference Krivov, Muff, Caflisch and Karplus2008; Rao & Caflisch, Reference Rao and Caflisch2004) with corresponding rates computed from the number of times a particular transition has occurred in the simulation. Alternatively, assuming that the structure is correlated with contact energy, one places edges between the nodes whose contact maps are sufficiently close (Hori et al. Reference Hori, Chikenji, Berry and Takada2009). From the network and estimated rate constants, a master equation can be constructed and studied analytically or solved numerically, using a method such as Gillespie simulations (Berezhkovskii et al. Reference Berezhkovskii, Hummer and Szabo2009; Buchete & Hummer, Reference Buchete and Hummer2008; Chodera et al. Reference Chodera, Singhal, Pande, Dill and Swope2007; Hinrichs & Pande, Reference Hinrichs and Pande2007; Hori et al. Reference Hori, Chikenji, Berry and Takada2009).

The master equation approach was developed in the context of protein folding, where the phase space is much larger compared with the size of the functional landscape. Sampling all of it with explicit solvent simulations is beyond computational reach for most of the current supercomputers. Therefore, work using this approach has to date mainly employed force fields with implicit solvent. However, with current state-of-the-art and future computational resources, it should be possible to accurately sample the protein's functional landscapes using explicit solvent.

Experiments by Frauenfelder and Ansari suggested that the basins of the native state ensemble contain sub-basins, which themselves may be subsequently divided further into substates, recursively, thus forming a hierarchical organization of the protein energy landscape (Frauenfelder et al. Reference Frauenfelder, Alberding, Ansari, Braunstein, Cowen, Hong, Iben, Johnson and Luck1990; Henzler-Wildman & Kern, Reference Henzler-Wildman and Kern2007). Based on this idea and somewhat related to the network approach above, the protein's landscape has been projected into a hierarchical tree with the help of a technique based mainly on PCA (Materese et al. Reference Materese, Goldmon and Papoian2008). This requires a long, straightforward sampling of the landscape in a simulation. PCA is a technique for analyzing trajectories based on coordinate transformation, where the new coordinates are the eigenvectors of the covariance matrix of the old coordinates. Thus, the new collective coordinates (called principal components) behave as independent random variables. They can be ranked by variance, so that the first principal component (PC) describes the largest fluctuations in the trajectory. In this way, most of the protein motions, particularly, large-scale motions (note that functional motions are usually large in scale), are well represented by only the first few PCs. When the simulation trajectory is projected into these PCs, multipeak distributions are observed, with peaks corresponding to different conformational states (García, Reference García1992). On the other hand, the vast majority of PCs are characterized by a narrow unimodal Gaussian-like distribution, corresponding to vibrational fluctuations (García, Reference García1992). A PC analysis of the 70-residue protein eglin c was recently carried out using the CHARMM force field with explicit solvent (Materese et al. Reference Materese, Goldmon and Papoian2008). The first two PCs, PC1 and PC2, contain the largest conformational fluctuations and turned out to be multipeaked. Several distinct regions of the (PC1, PC2) plot are highly populated, or, represent basins of the free energy surface (Fig. 7).

Fig. 7.

Tree-like hierarchy of the functional landscape of eglin c is revealed by PCA. In the higher index PCs, the basins show finer splitting, like a basin on a (PC2, PC3) plot that splits into two sub-basins on a (PC3, PC4) plot. Five hierarchical tiers were identified. The number of tiers up to the first common parent basin can serve as ultrametric ‘distance’ between any two sub-basins on the same tier. This ultrametric can be used as a basis to generalize rate constants in kinetic schemes of landscape mapping from partial sampling of the landscape (Materese et al. Reference Materese, Goldmon and Papoian2008). We thank Christopher K. Materese for preparing the figure.

Furthermore, zooming into a particular basin allows exploration of finer splittings: plotting the population using PCs of higher number as coordinates reveals clearly identifiable sub-basins. Each of them may be scrutinized further in the PCs of even higher number and so on if there is sufficient sampling of the next-level sub-basin from the MD trajectory. Hence, the novelty of this approach was not only to project the protein dynamics into the first two PCs, which is common, but also follow basin splitting in PCs of higher dimensions (Materese et al. Reference Materese, Goldmon and Papoian2008). For eglin c, five hierarchical tiers were found from a 330 ns simulation. Figure 7 shows a tree-like organization of conformational states, providing microscopic detail to earlier conceptual arguments on the hierarchical organization of the protein native state (Frauenfelder et al. Reference Frauenfelder, Alberding, Ansari, Braunstein, Cowen, Hong, Iben, Johnson and Luck1990; Henzler-Wildman & Kern, Reference Henzler-Wildman and Kern2007).

The tree plays a role similar to that of networks discussed above, that is, shows the connections between the conformational states of the same tier. The tree defines an ultrametric (An ultrametric is a metric that satisfies the following strengthened version of the triangle inequality. Metric is a rule that defines a ‘distance’ between any two elements of the set. If it is an ultrametric, then among three pairwise distances between any three elements, two are equal and the third one is smaller. An ultrametric in a tree can be defined as the number of levels upward to the first common ancestral branch.) among its same-level branches, and assuming that the ultrametric reflects actual physical properties of the energy landscape, such as barrier heights, one can use it to assign the transition rates between the conformational states on a particular tier: the closer the states by the tree ultrametric, the higher the transition rate between them. In the next approximation, typical barrier heights or even rates can be calculated for each possible ultrametric distance (there is only as many values as there are levels in the trees) by performing separate simulations, several for each distance. This would be a way for the reconstruction of the energy landscape ‘from the bottom’.

This discussion of the populated basins in terms of PCs suggests that PCs may be good coordinates for yet another representation of the protein energy landscape – construction of free energy surfaces (García & Sanbonmatsu, Reference García and Sanbonmatsu2001; Maisuradze et al. Reference Maisuradze, Liwo and Scheraga2009b). They have been used for that purpose in combination with the kinetic network approach, for protein G and SH3 domain (Hori et al. Reference Hori, Chikenji, Berry and Takada2009). In that study, PCs were based not on Cartesian coordinates of the atoms or dihedral angles, but on the contact order map of the protein, taking advantage of the funneled shape of the energy landscape. Along the same lines, Q – the fraction of contacts shared with the native state – may also be used. That approach has been used for 20-residue peptide Trp-cage (Wu et al. Reference Wu, Zhuravlev and Papoian2008). The authors took two native-like structures from an NMR study and defined two Qs, with respect to either of them, using these two coordinates to build a 2D free energy surface, unfurling the native basin at high resolution (Wu et al. Reference Wu, Zhuravlev and Papoian2008), which is elaborated below.

4.2 Applications of the landscape paradigm to protein dynamics and function

When a landscape representation is available, detailed investigation may shed light on the nature of protein dynamics and function, such as mechanisms of enzymatic catalysis or molecular motor stepping. A free energy surface in two dimensions is a potent representation of the underlying energy landscape, although building it requires significant computational efforts (Pisliakov et al. Reference Pisliakov, Cao, Kamerlin and Warshel2009; Wu et al. Reference Wu, Zhuravlev and Papoian2008). As a first step, one has to choose good order parameters as two collective coordinates – that is, the remaining transverse degrees of freedom should be much less important for the questions asked, and, preferably, quickly equilibrate. These are strong requirements, and a choice of coordinates, therefore, is a significant approximation. Would a pair of coordinates be enough for an adequate description? How many coordinates does one need? As an example, Materese et al. have shown through sub-microsecond simulations that at least in one case, eglin c, only on the order of 5 out of about 300 conformational coordinates are required for a comprehensive description of the non-trivial part of the native protein dynamics (Materese et al. Reference Materese, Goldmon and Papoian2008). However, not all features of the protein's conformational dynamics might be of importance for a given problem; hence, lower-dimensional surfaces still play a significant role.

PCA relies on the actual dynamics of a protein system during some period of observation to find the coordinates that optimally represent motions. If there is no initial knowledge about the system dynamics or energy landscape, the choice of collective coordinates must be based on physical plausibility. Q is correlated with folding, but in the bottom of the funnel the motions transverse to Q may also be very important because of the rugged topography of the functional landscape. Interestingly, the choice of two Qs, with respect to two different (but similar) basins in the native state, provides both correlation with the funneled shape of the landscape, as a Q ₁=Q ₂ diagonal of the (Q ₁=0 … 1, Q ₂=0 … 1) square, and description of the motions transverse to folding, which are transverse to the above-mentioned diagonal (Wu et al. Reference Wu, Zhuravlev and Papoian2008). The diagonal direction allows describing native excitations that occur through partial unfolding. The transverse to diagonal motions are resolved non-uniformly throughout the square – they are highly resolved as the native state is approached in the upper-right corner (see Fig. 8) and completely unresolved near the origin (Wu et al. Reference Wu, Zhuravlev and Papoian2008). After computing Free energy surface (FES) with the chosen coordinates, Langevin dynamics may be used to run trajectories on this surface, which requires many orders of magnitude less computational effort than a direct simulation of protein dynamics. Therefore, it becomes feasible to run tens of thousands of trajectories, producing various first passage time moments and distributions.

Fig. 8.

A free energy surface is constructed for Trp-cage through 2D umbrella sampling. The two coordinates are the fractions of contacts shared with either of two similar structures labeled as ‘1’ and ‘37’. Magenta and white on the upper panel denote two overdamped Brownian dynamics trajectories. The upper right panel shows a magnified region of the surface, revealing the sensitivity of the trajectory to fine features of the surface, such as shallow traps and low ridges. These features broaden the distribution of first passage times from ‘37’ to ‘1’, but there are no conspicuous deep traps that would engender classification. To get reference for the first passage time distribution at room temperature, two more First-passage time (FPT) distributions were calculated at other temperatures with their coefficients of variation–square-root-of-variance-to-mean ratios (lower left panel). C _v for the room temperature is near 1, which is considered a boundary between ‘narrow’ and ‘broad’ distributions. The dependence of the mean first passage time on temperature also changes around room temperature from Arrhenius (corresponding to activated dynamics) to reverse proportionality (corresponding to confined diffusion) (lower right panel). These facts lead to conclusion about the dynamical behavior of Trp-cage at room temperature being borderline between activated and diffusional (Wu et al. Reference Wu, Zhuravlev and Papoian2008).

One important question about native basin energy landscapes concerns the nature of conformational dynamics occurring on it. In particular, dynamics are called activated (Debenedetti & Stillinger, Reference Debenedetti and Stillinger2001) if the molecule spends most of the time in the landscape minima, rarely hopping from a minimum to a minimum, with the average waiting time given by the Arrhenius law. This is a typical situation for a liquid that is slightly supercooled. When the relaxation times become longer than the observation window, the dynamics are considered glassy. On the contrary, if the system spends most of the time cruising from one saddle point to another on the landscape, the dynamics are diffusive, similar to those of a normal liquid (La Nave et al. Reference La Nave, Scala, Starr, Sciortino and Stanley2000; Schulz et al. Reference Schulz, Krishnan, Daidone and Smith2009).

To ascertain which of these dynamical regimes is relevant for a small protein Trp-cage, Wu et al. computed a 2D free energy surface using two structure-based order parameters, and ran on it hundreds of thousands of Brownian dynamics trajectories (Wu et al. Reference Wu, Zhuravlev and Papoian2008). In particular, different sets of simulations were computed at various temperatures. By analyzing the temperature dependence of the mean first passage times to transition from one native-like state to another as well as the shapes of the first passage times distribution, it was concluded that Trp-cage room temperature dynamics are borderline between diffusive and activated dynamics (Wu et al. Reference Wu, Zhuravlev and Papoian2008). To further refine these types of simulations, the dependence of the friction constant on the surface position should be taken into account (Socci et al. Reference Socci, Onuchic and Wolynes1996; Wang et al. Reference Wang, Zhang, Lu and Wang2005). The latter inhomogeneity arises from the ruggedness created by the degrees of freedom that are transverse to Q ₁ and Q ₂. However, these types of simulations are most successful when surface coordinates are appropriately chosen to describe a particular function or process. For instance, Heath et al. have elucidated the mechanisms of misfolding in a mutant of a protein S6 with the help of a free energy surface in terms of Q (fraction of native contacts) and the fraction of non-native contacts (Heath et al. Reference Heath, Kavraki and Clementi2007).

In a different example, where the need for two collective coordinates arose naturally, it was recently suggested from interpretation of experiments that allosteric and catalytic motions are coupled in adenylate kinase (Henzler-Wildman et al. Reference Henzler-Wildman, Lei, Thai, Kerns, Karplus and Kern2007a, Reference Henzler-Wildman, Thai, Lei, Ott, Wolf-Watz, Fenn, Pozharski, Wilson, Petsko, Karplus, Hubner and Kernb; Kern et al. Reference Kern, Eisenmesser and Wolf-Watz2005). To address this question in detail, Pisliakov et al. performed multiscale simulations (Pisliakov et al. Reference Pisliakov, Cao, Kamerlin and Warshel2009), building a 2D free energy surface in coordinates relevant to the question: a chemical coordinate (corresponding to the chemical reaction being catalyzed) and a conformational coordinate (corresponding to the conformational motion of the enzyme). Pisliakov et al. performed MD simulations on two scales: all-atom and coarse-grained, changing the parameters of the coarse-grained model such that it reproduces the autocorrelation functions computed from atomistic dynamics. Then the analytical 2D surface was built, and the parameters were tuned in a similar fashion to achieve agreement with both of the more detailed models and also two rate constants obtained in FRET and NMR studies (Henzler-Wildman et al. Reference Henzler-Wildman, Thai, Lei, Ott, Wolf-Watz, Fenn, Pozharski, Wilson, Petsko, Karplus, Hubner and Kern2007b). Studying the dependence of the first passage times from the reactant state to the product state at various barrier heights and friction constants, Pisliakov et al. have concluded that their simulations indicate no direct channeling of conformational motion excitations into chemical transformation (Pisliakov et al. Reference Pisliakov, Cao, Kamerlin and Warshel2009). This was explained by a timescale separation between picosecond and nanosecond solvent motions and coherent protein modes, on the one hand, and millisecond chemical catalysis step, which is long enough for the coherent protein motions to completely dissipate (Pisliakov et al. Reference Pisliakov, Cao, Kamerlin and Warshel2009). These types of questions may be explored in the future with the help of ab initio MD simulations, which treat protein motions and chemical rearrangements on an equal footing (Papoian et al. Reference Papoian, Ulander and Wolynes2003b).

When function is strongly coupled to the folding landscape, as may be the case with intrinsically disordered proteins (Papoian, Reference Papoian2008), coordinates related to folding have to be used for building the FES. Using Qs, Turjanski et al. (Reference Turjanski, Gutkind, Best and Hummer2008) have constructed an FES for a natively unstructured transcription factor, and found that binding of the ligand precedes folding of the transcription factor, indicating a coupling between binding and folding in accordance with NMR experimental data (Sugase et al. Reference Sugase, Dyson and Wright2007). In a computational follow-up of an interesting recent finding of catalytic activity in a disordered enzyme (Pervushin et al. Reference Pervushin, Vamvaca, Vogeli and Hilvert2007), radius of gyration and contact order were used for the 2D free energy surface (Roca et al. Reference Roca, Messer, Hilvert and Warshel2008). Subsequent analysis indicated that upon substrate binding the molten globule-like phase partially folds (Roca et al. Reference Roca, Messer, Hilvert and Warshel2008). Coupling between binding and folding is a known and frequent phenomenon (Papoian & Wolynes, Reference Papoian and Wolynes2003). Evolution may have produced natively disordered enzymes to provide finer temporal control over enzymatic activity (as unfolded chains are faster degraded by proteosomes) and acceleration of the binding kinetics through a ‘fly-casting’ mechanism (Hyeon & Onuchic, Reference Hyeon and Onuchic2007; Shoemaker et al. Reference Shoemaker, Portman and Wolynes2000) (as an unfolded chain occupies a larger volume) (Papoian, Reference Papoian2008). One of the possible realizations of the coupling between binding and folding is a frustration localization near the binding site. Upon binding, the frustration disappears and induces a conformational change, resulting in folding (Ferreiro et al. Reference Ferreiro, Hegler, Komives and Wolynes2007). Papoian and Wolynes developed a theory of the coupling between binding and folding energy landscapes (Papoian & Wolynes, Reference Papoian and Wolynes2003), which can be used to study binding and folding thermodynamics and kinetics (Levy et al. Reference Levy, Papoian, Onuchic and Wolynes2004; Marcovitz & Levy, Reference Marcovitz and Levy2009; Wang et al. Reference Wang, Zhang, Lu and Wang2006). Currently, computational power allows construction of the 2D FES even with explicit solvent all-atom force fields, as was applied to study metal-coupled folding of a zinc-finger motif (Li et al. Reference Li, Zhang, Wang and Wang2008).

With the help of free energy surfaces, Hyeon and Onuchic have elucidated the mechanism of chemical energy transduction into mechanical motion in kinesins (Hyeon & Onuchic, Reference Hyeon and Onuchic2007). Kinesin is a molecular motor that walks in 8 nm directional steps along microtubules consisting of 13 protofilaments. Kinesin consists of two heads, connected by a neck region, that bind to specific sites on the protofilaments. A conformational switch of the neck, called a power stroke, precedes binding of the tethered (not bound) head to the microtubule, which, in turn, occurs after a non-directional diffusive search for a binding site. In this work, several different free energy surfaces were computed for two different purposes (Hyeon & Onuchic, Reference Hyeon and Onuchic2007). First, the authors considered the question of which microtubule-binding site is preferred by a motor head. To address it, Hyeon and Onuchic constructed an FES corresponding to a potential of mean force for the tethered head as a function of its three spatial coordinates. Two such surfaces were computed for both conformational states of the neck region using hybrid Monte Carlo–MD simulations with a coarse-grained force field (Hyeon & Onuchic, Reference Hyeon and Onuchic2007). Examination of the surface cross-sections allowed the researchers to identify several basins of attraction corresponding to specific binding sites on the same or different protofilaments in addition to a separate broad entropic basin. The most favorable binding site turned out to be on the same protofilament that hosts the other motor head. Hyeon and Onuchic then combined the two surfaces (for two different conformational states of the neck region) in a single time-dependent surface that would imitate the conformational switching of the neck region during the power stroke. In this combined potential, they ran Brownian dynamics treating the tethered head as a quasi-particle. Depending on the timescale of potential switching (the power stroke), the trajectory may be trapped in the intermediate state, corresponding to the head bound to the nearest binding site (the one on the neighboring protofilament, see Fig. 9). Hyeon and Onuchic also built free energy surface in terms of Q _p and Q _int, which are the fraction of native contacts in the microtubule-binding motif of kinesin and the fraction of interfacial native contacts between kinesin and microtubule. Examining these surfaces showed that cracking – disruption of some kinesin native contacts – occurs prior to binding. The cracking increases the binding rate because of enhanced flexibility of the partially unfolded structure, through the fly-casting mechanism. Interestingly, the cracking happens right before the moment of binding.

Fig. 9.

A microtubule consists of 13 protofilaments made of tubulin heterodimers. Kinesin motors walk on microtubules directionally in 8 nm steps. Hyeon and Onuchic have shown that when one motor head is bound the other one binds to the active site on the same protofilament with a possible intermediate state on the neighboring binding site. (Hyeon & Onuchic, Reference Hyeon and Onuchic2007).

Motivated by the allostery problem, the idea of cracking may also be used to extend the Gō model to a protein that has two well-defined conformational states with a transition between them. This method of constructing protein energy landscapes was used to study glutamine-binding protein (GBP), S100A6 (structural analogue of calmodulin), dihydrofolate reductase and HIV-1 protease (Okazaki et al. Reference Okazaki, Koga, Takada, Onuchic and Wolynes2006). Let us assume that there is a Gō-like force field V ₁(R, R ₁), where R represents the coordinates of an arbitrary conformation and R ₁ are the coordinates in the native conformation. If there is another native basin, the corresponding Gō model will be given by a potential V ₂(R, R ₂). Next, one wants to construct a single potential where both basins are present, though, perhaps with a difference in stability (for example, the second basin is shifted by ΔV). To merge the two surfaces together, techniques from the Marcus theory of electron transfer (Marcus & Sutin, Reference Marcus and Sutin1985) may be borrowed, and used to ‘sew’ the two Gō-potentials together. Introducing a coupling constant Δ, Okazaki et al. (Reference Okazaki, Koga, Takada, Onuchic and Wolynes2006) obtained a smooth double-basin potential V _MB given by

the equation for eigenvalues. The eigenvector identifies the basin in which the protein is residing and therefore may be used to construct a reaction coordinate. Plotting free energy of GBP as a function of this reaction coordinate in addition to its components – the average energy and conformational entropy – led the authors to observe a sudden increase in conformational entropy in the middle region of the transition profile. This is the hallmark of cracking – partial unfolding and breaking some of the native contacts (Hyeon & Onuchic, Reference Hyeon and Onuchic2007; Hyeon et al. Reference Hyeon, Jennings, Adams and Onuchic2009; Miyashita et al. Reference Miyashita, Onuchic and Wolynes2003, Reference Miyashita, Wolynes and Onuchic2005; Okazaki et al. Reference Okazaki, Koga, Takada, Onuchic and Wolynes2006; Whitford et al. Reference Whitford, Miyashita, Levy and Onuchic2007, Reference Whitford, Gosavi and Onuchic2008a, Reference Whitford, Onuchic and Wolynesb). Again, it turned out to be useful to plot free energy surfaces in terms of Q ₁ and Q ₂, that is, the fractions of contacts shared with either one of the native states, respectively. For GBP this surface showed that there is a straight valley between the two minima corresponding to the two states, which means that breaking of contacts specific to the initial basin and formation of contacts specific to the final basin occur simultaneously.

The allosteric transition of S100A6 is more sequential: first, contacts specific to the first basin break and then new contacts are formed, leading to greater unfolding. Different coordinates for the surfaces, Q _common and Q ₁–Q ₂, where Q _common is a fraction of the contacts common to both basins that are present in the conformation, provided another perspective on the transition dynamics. Authors plotted free energy surfaces at the folding temperature and then below the folding temperature and observed the appearance of an intermediate state for S100A6 at the lower temperature (Okazaki et al. Reference Okazaki, Koga, Takada, Onuchic and Wolynes2006). Another method of combining several Gō potentials, based on summation of the partition functions, rather than secular equation, has also been reported (Best et al. Reference Best, Chen and Hummer2005).

A different way to conjugate two structure-based (Gō-like) potentials has been used with adenylate kinase (Whitford et al. Reference Whitford, Gosavi and Onuchic2008a). The enzyme consists of three domains: CORE, NMP and LID. The conformational change from open to closed structure consists of motions of the NMP and LID domains. In this study, a potential based on the open structure was modified by the addition of a bias favoring contacts present in the closed structure. A 2D free energy surface was constructed through MD simulations in terms of distances between CORE and NMP centers-of-mass and between CORE and LID centers-of-mass (luckily, adenylate kinase has these natural reaction coordinates). Whitford et al. also constructed transition paths for opening and closing using the adiabatic (non-linear) normal mode technique discussed in the first section of the review (Miyashita et al. Reference Miyashita, Onuchic and Wolynes2003). Interestingly, the pathways largely follow through the valleys of the landscape, especially, the closing transition.

Yet another possible set of coordinates to choose for plotting free energy surfaces may be derived through the Scalable ISOMAP technique (Das et al. Reference Das, Moll, Stamati, Kavraki and Clementi2006). From a practical viewpoint, it appears similar to PCA, that is, it reduces the dimensionality of the dynamics, but does so in a non-linear way. Using surfaces with these coordinates, Haspel et al. identified contacts that shape the native basins of C3d/Efb-C and two of its mutants (Nurit et al. Reference Nurit, Geisbrecht, Lambris and Kavraki2010).

4.3 Role of water in shaping the landscape

Many of the works in this section, especially those dealing with large proteins, relied on implicit solvent Hamiltonians (Chen & Brooks, Reference Chen and Brooks2007; Pincus et al. Reference Pincus, Cho, Hyeon and Thirumalai2008). However, the role of water in sculpting the functional landscape is extremely important, as elucidated by Materese et al. (Reference Materese, Goldmon and Papoian2008) through examining water-mediated contacts, which may cause splitting of native basins into sub-basins. The presence of water-mediated contacts is frequent in the native state (Zong et al. Reference Zong, Papoian, Ulander and Wolynes2006; Papoian et al. Reference Papoian, Ulander, Eastwood, Luthey-Schulten and Wolynes2004). The role of water in protein folding has also been appreciated for a long time (Cheung et al. Reference Cheung, García and Onuchic2002; Head-Gordon & Brown, Reference Head-Gordon and Brown2003; Kauzmann, Reference Kauzmann1959; Kaya & Chan, Reference Kaya and Chan2003; Lucent et al. Reference Lucent, England and Pande2009; van der Vaart et al. Reference van der Vaart, Bursulaya, Brooks and Merz2000), where the hydrophobic effect is considered to be the primary driving force for folding (Dyson et al. Reference Dyson, Wright and Scheraga2006). However, water-mediated hydrophilic interactions were also found to play an important role (Papoian et al. Reference Papoian, Ulander, Eastwood, Luthey-Schulten and Wolynes2004). Dynamic coupling between water and proteins usually occurs on fast timescales. When a protein is effectively caged by solvent, the motions of the former become slaved to those of the latter (Austin et al. Reference Austin, Beeson, Eisenstein, Frauenfelder and Gunsalus1975; Frauenfelder et al. Reference Frauenfelder, Fenimore and Young2007). Such coupled collective modes can be probed by THz spectroscopy (Ebbinghaus et al. Reference Ebbinghaus, Kim, Heyden, Yu, Gruebele, Leitner and Havenith2008; He et al. Reference He, Ku, Knab, Chen and Markelz2008).

On the other hand, functional motions of a protein chain are slow, so there is a timescale gap between the protein and solvent motions. Therefore, in most of the problems discussed in this review, the picture of a chain moving in a thermally averaged field of solvent motions is rather appropriate. In this way, water smoothes the energy landscape (Zhuravlev et al. Reference Zhuravlev, Materese and Papoian2009), which, in turn, allows us to talk about effective interaction energies, where solvent degrees of freedom are integrated out. Although implicit solvent models are physically justified by these considerations and are frequently used, it is desirable to account for the role of water as accurately as possible, for instance, by including water-mediated contacts through a second well in potentials between polar groups, or through the introduction of non-pairwise potentials that are dependent on residues' water exposure (Papoian et al. Reference Papoian, Ulander, Eastwood, Luthey-Schulten and Wolynes2004; Zong et al. Reference Zong, Papoian, Ulander and Wolynes2006). The resulting representation of the long-range hydrophilic contacts is more rigorous and improves structure prediction (Papoian et al. Reference Papoian, Ulander, Eastwood, Luthey-Schulten and Wolynes2004; Zong et al. Reference Zong, Papoian, Ulander and Wolynes2006). Alternatively, explicit solvent simulations may be used, which are computationally much more demanding.