Gō model revisited.

This review discusses Gō models broadly used in biomolecular simulations. I start with a brief description of the original lattice model study by Nobuhiro Gō. Then, the theory of protein folding behind Gō model, free energy approaches, and off-lattice Gō models are reviewed. I also mention a stringent test for the assumption in Gō models given from all-atom molecular dynamics simulations. Subsequently, I move to application of Gō models to protein dynamical functions. Various extension of Gō models is also reviewed. Finally, some publicly available tools to use Gō models are listed. 2019 © The Biophysical Society of Japan.

Since the very first paper by Taketomi, Ueda, and Gō in 1975, the so-called Gō model has been broadly used in various biomolecular simulations. A brief history, a stringent test, various extensions, and broad range of applications are reviewed.

Since the very first paper by Taketomi, Ueda, and Gō in 1975 [1], the so-called Gō model has long been used to broad range of biomolecular simulations, primarily for, but not limited to, protein folding studies. In this review article, I start with a brief description of the original work by Gō, which is followed by subsequent developments, more recent application of Gō model, and discussion of future directions.

While not many, in the last 20 years, the Gō model has been reviewed in a few articles [2-4], which would complement this review.

Lattice model by Gō

Teketomi, Ueda, and Gō, for the first time, introduced a lattice model to protein folding study [1]. In their model, each monomer in “protein” is placed on a lattice point in two dimension and is connected by bonds that have the unit length of the lattice (Fig. 1A). Starting from random configurations, the Metropolis Monte Carlo simulation was used to fold this model protein. Monomers that are separate by the unit length and not connected by a bond have non-local contact interactions. For the contact interactions, they considered three cases, a weak limit of specificity, an intermediate specificity, and a strong limit of specificity. It is in the third class that turned out to be called Gō-model later; two beads have negative contact interaction energy only when these pairs are in contact at the pre-defined native structure. Otherwise, the two beads have no contact energy. The Monte Carlo simulations resulted in complete folding to the native structure only for the strong limit of specificity, but not for the other two cases. A subsequent paper by Gō and Taketomi introduced the specificity to the local potential, a negative energy only when the bond angle is the same as that in the predefined native one, studying the balance between local and non-local specific interactions [5]. In these and other studies, Gō and his collaborators pursued statistical physics of the strong limit of specificity, finding many qualitatively consistent results with experiments, as summarized in the seminal review [6]. In their lattice model, the protein folds to the native configuration below the melting temperature, above which the protein unfolds. The transition is fairly cooperative, showing a clear peak in the heat capacity. The folding occurs faster when the local preference energy is larger, relative to the non-local contact energy, whereas the folding transition is more cooperative when the non-local contact contribution is larger. All these are very robust physical properties well supported by biophysical experiments and shared by more accurate protein simulations performed to date.

Figure 1

Various studies with Gō models. A) One of the original two-dimensional lattice models used by Gō. The figure taken from [5]. B) A stringent test of the assumption in Gō models. For contact i and j, the log ratio of lifetime in the transition-path tTP to lifetime in the unfolded state tU is plotted by color in upper-left triangle, while the native contact map is depicted in the lower-right triangle. Results for three proteins are drawn here. The figure taken from [35]. C) Comparison of the root mean square fluctuation in the native basin for a test protein CheY. Results from all-atom MD (black circles), an elastic network model (GNM) (red curve), and a Gō model (green curve) are compared. The figure taken from [30]. D) A schematic plot for the multiple-basin Gō model that is based on two single Gō models. The picture taken from [41].

About 15 years later from the first paper by Gō, the lattice model of protein folding got popularity, notably by the studies of Dill, Shakhnovich, and Onuchic and Wolynes [7-10]. Importantly, in this second generation of the lattice protein model, monomer contact energies were not biased to their pre-defined native structure, but are purely sequence-based. So, these models are not a sort of Gō model. But, sometimes, the Gō model was mentioned as a control, which served as the limit of strong specificity. It was in this context that people started to call this type of models as “Gō model” [7,8]. Retrospectively, due to its strong geometrical constraint, the sequence-based, i.e., non-Gō type lattice models may have pronounced ruggedness, which might have overestimated the ruggedness of the energy landscape, comparing with nowadays-available all-atom molecular dynamics (MD) simulations for folding.

Theory behind

Through series of lattice model simulations and lessons from many X-ray crystallographic structures, Nobuhiro Gō proposed the consistency principle in the 1983 review [6]. In their native structures, proteins individually take the optimal local interactions, akin to the secondary-structure, and non-local interactions, i.e., the tertiary structure. This dual optimality was attained via evolutionary selections of the amino acid sequence. The review also mentioned about the kinetic aspect of the consistency principle; the balance between the local and non-local specificity, which is necessary for the well-behaved cooperative folding-unfolding transition. The consistency principle serves as a theoretical basis to use the Gō model. The consistency principle, and thus the Gō model, were considered as an ideal limit, whereas real proteins must have other restraints, such as functional restraints.

Some years later, extending Gō’s perspective, Wolynes and his coworkers formulated a statistical physical theory of protein folding and elucidated global picture of energy landscape of proteins [11,12]. The theory clarified that, for a given temperature, the energy bias to the native structure has to be large enough, relative to the average ruggedness of the energy surface (frustration) for the protein being able to fold to the native structure avoiding glassy slow dynamics. Here, the energy bias is linked to the thermodynamic stability T < Tm, where T is the physiological temperature and Tm is the melting temperature. The average ruggedness is correlated with the characteristic glass-transition temperature Tg for the onset of slow dynamics; thus, as the kinetic condition, T > Tg is to be satisfied. Together, Tg < T < Tm gives the foldability condition. For fast folding of proteins, the frustration at the native basin and on the route to it has to be small enough, which is termed the principle of minimum frustration.

While the consistency principle of Gō and the principle of minimum frustration apparently have overlap, they also differ in some respects. First, the consistency principle is a concept stated primarily by words, while the principle of minimum frustration gives mathematical expressions. Secondly, the consistency principle primarily states on the structural aspect at the native state and thus on the thermodynamic stability, in addition to the balance between the local and non-local specificities for experimentally observed cooperative folding-unfolding transitions. The Wolynes theory defines the foldability as the combination of the stability and kinetics; T < Tm from the stability and T > Tg from the kinetics.

Subsequent computer simulations by Onuchic and others pointed out importance of the network of folding pathways [13]; the folding is fast and efficient when there are multiple parallel pathways that are linked to the native state, while it is slow when there is a bottleneck where only few routes exist. This ends up with the concept of protein folding funnel; fast folding proteins have funnel-like energy landscapes [14]. The solvent-averaged energy of a fast-folding protein, on average, decreases as the protein approaches to the native state, which correspond to the bottom of the funnel. If the effective energy monotonically decreases, how does the folding free energy barrier arise? As the folding reaction proceeds, the number of available conformations, and thus the conformational entropy decreases, which opposes the folding. It is thus the energy-entropy compensation that gives rise to the free energy barrier along the folding reaction coordinate. Conformational bottleneck corresponds to a large conformational entropy loss, which increases the free energy at the bottleneck.

Gō model can be considered as a concise mean to realize the consistency principle, or the perfect funnel landscape. Assumed in Gō model is that only natively-contacting pairs (native contacts) have crucial contribution to the folding, but other pair interaction (non-native contacts) can be ignored.

Free energy function type Gō model

In the same spirit of Gō model, Wako and Saito proposed an Ising-like model to describe protein folding processes [15,16], which was, many-years later, reinvented and applied to many proteins by Munoz and Eaton [17]. In the model, each residue takes two state, n (native) and u (unfolded) so that a microstate of a protein is described as a sequence of n and u, such as uuuunnnnuu for a 10-residue case. The global native state N and unfolded state U correspond to the states that all the residues take n and u microstates, respectively. Only when both of natively interacting residues as well as all the residues in between take the n states, they get a certain stability, which basically is the same assumption as the Gō model. Non-native contacts are completely ignored. As residues change from the u state to the n state, the chain loses its conformational entropy. Thus, as the protein folds, it gains energetic stability and loses its conformational entropy. Since the folding transition is modeled as the energy-entropy compensation process, the optimal folding pathways are determined so as to minimize the entropy loss for a given energy gain. Under some assumptions, free energy of microstates can be analytically calculated by a transfer-matrix approach, which is a clear advantage of this type of modeling, compared to the lattice Gō model. This simple modeling turned out to have marked predictive power of folding transition ensembles and folding rates, when compared with experimental phi-value analysis and folding rates. Related free energy function approaches were developed by a few other studies, as well [2,18-20]. Conversely, this success strongly supports the consistency principle and the perfect funnel view. More recently, the Wako-Saito-Munoz-Eaton model and its extension have been broadly used in diverse folding studies [21,22].

Off-lattice Cα Gō model

While the lattice model was powerful to reveal conceptual aspects in protein folding, its highly restricted geometry precludes direct application to real protein structures. Especially, low-frequency and collective fluctuations around the native structures cannot easily be represented by lattice models. In this respect, another class of minimal protein models, called the elastic-network-model, was invented first by Tirion in 1996 [23], and further developed by Jernigan, Bahar, and others [24]. The elastic network model consists of elastic bonds between two monomers that are spatially close to each other at the native structures, with the natural lengths of the elastic bond being the lengths at the native structure. Albeit its ultra-simplicity, the elastic network model was shown to represent low-frequency fluctuations of proteins in the native state surprisingly well. As in the case of Gō model, the elastic network model directly uses and is biased towards the native structure by construction. However, because all the interactions are elastic, the elastic network model cannot approximate the unfolded state at all, and thus is not considered as a Gō model.

The lattice Gō model is a concise realization of protein folding, but not good for native fluctuation dynamics. On the other hand, the elastic network model is a simple and minimal model to approximate low-frequency native fluctuations, but does not take into consideration of unfolding. Soon after Tirion’s work, Clementi, Nymeyer, and Onuchic proposed an off-lattice Gō model that resulted in taking advantage of the two minimal models; the model represents a protein as a chain of Cα atoms of every amino acids and that has both local angle and non-local contact potentials biased towards the native structure. In the low temperature limit, the model converges to the elastic network model, whereas the protein unfolds cooperatively at a higher temperature. They applied this off-lattice Cα Gō model to small fast-folding proteins, directly comparing the folding pathways with experiments [25]. Several subsequent works together showed its predictive power for folding reaction mechanisms [26,27].

As is clear by now, “Gō model” does not mean a single model, but represents a class of models that share the concept in the original work by Gō. Although no one clearly defined it, in practice, Gō models share the two concepts. 1) They represent folding-unfolding transitions and 2) take into accounts pairwise contact energies only for the pairs that are in contact in the native structure. The Gō model is often called “structure-based model” as well since the energy function is explicitly dependent on the native structure. While the two names are used as synonym, they can be slightly different; e.g., the elastic network model of Tirion is clearly a structure-based model, but is not a sort of Gō model.

The first version of off-lattice Cα Gō model by Clementi, C., et al. is a minimal and concise model so that there is room to add some more details. Karanicolas and Brooks added sequence-dependent local terms, chemically-motivated energy scales for the contact energies, and a desolvation potential in non-local contacts, improving predictive power of folding transitions and pathways [28]. A similar desolvation energy was addressed by Kaya and Chan, as well [29]. Multiscale algorithms were utilized to reflect more sequence-and structure-specific interaction at atomic level within Cα Gō model [30,31]. Furthermore, the chemical denaturant effect was incorporated into an off-lattice Gō model [32].

Testing Gō model assumption via all-atom molecular simulations

A special-purpose supercomputer Anton and the associated MD software Desmond were used to simulate long-time folding dynamics of small proteins with an all-atom force-field and explicit water solvent, which successfully realized reversible folding and unfolding for more than 10 proteins [33,34]. They observed in the trajectories that, whereas nonnative secondary structures do form in the unfolded state ensemble for some proteins, they are transient and disappear either before initiation of folding events or early during the folding. This implies that these non-native contacts do not contribute much to the folding mechanisms, supporting the assumption behind Gō models.

Soon after this all-atom folding simulations, these trajectories were further analyzed from the perspective of folding mechanisms by Best, Hummer, and Eaton [35]. Probably, this analysis provides the most convincing data to date that supports the native-contact centric view and thus the use of Gō model. For each residue pair, they compared the lifetime of that contact in the transition-path and the lifetime of the same contact in the non-native ensemble in the trajectories (the transition-path stands for a fragment of MD trajectory that departs from the unfolded basin and reaches at the native state) (Fig. 1B). If the former lifetime is much longer than the latter for a pair, this indicates the pair has significant role for folding. Quantifying this ratio, they found that the high scores are located only in the native contact pairs and their neighbors for all but one protein analyzed. Interestingly, one exception was a designed protein, of which sequence has not evolved naturally. For naturally-evolved proteins, not all the native contacts are equally resistant; the high scored regions correlate with the folding initiation sites. They also performed a Bayesian analysis finding similar tendency; native contacts are formed at higher probabilities in the transition-path than the non-native contacts. This analysis unambiguously shows that, in the successful folding and unfolding transitions, natively formed contacts play major roles, while non-native contacts do not contribute significantly; thus, this directly rationalizes the use of Gō models in predicting/revealing folding mechanisms for natural proteins. Furthermore, the Wako-Saito-Munoz-Eaton model was applied to the same set of proteins as those simulated with Anton, finding highly consistent results [36].

From folding to function with Gō model

It has been shown that off-lattice Gō models are useful not only for protein folding, but also for native-state dynamic simulations [30,37]. Somewhat surprisingly, comparing with the root-mean-square-fluctuations (RMSF) calculated by all-atom MD with explicit water, off-lattice Gō models agree better than the elastic network model, which indicates that even the native dynamics of proteins reflect some un-harmonic part of potentials that is linked to local unfolding (Fig. 1C). It was investigated that, relative to the elastic network model, local unfolding, or cracking, decreases the free energy barrier for conformational transition [38]. With a multiscale-calibrated version of Gō model, the correlation between the RMSF by the Gō model and that by all-atom MD with explicit water was as high as the correlation between the RMSF by all-atom MD with implicit water and that with explicit water [39].

Given that the native dynamics can be well captured by Gō models, we can go further simulating conformational transition between two or more distinct conformations relevant to allosteric regulations. When two distinct conformations, A and B, for a target protein are available from experiments, we can construct the two respective Gō models V(R|A), and V(R|B), each making funnel-like landscapes centered around A and B basins, respectively, where the vector R collectively represents all the particle coordinates. If we think of V(R) = Min (V(R|A), V(R|B)), this potential V(R) encodes two basins A and B in its energy landscape (Fig. 1D). To use it in MD simulations, we need to make it differentiable. There can be multiple means to do it. Best and Hummer proposed the so-called soft-min function, represented as a combination of logarithmic and exponential functions, describing conformational change of a protein [40]. Okazaki, K., et al. used the secular equation formalism for the purpose [41]. Notably, these multiple-Gō model formalisms create two basins in the global energy landscape, but do not make intermediate states between two conformations. As an opposite limit of cooperativity, one can put two local minima in every pair contact potentials, corresponding to the pair distances in conformations A and B [42,43]. While this formalism also encodes two basins A and B, it also allows to possess multiple intermediate states between two endpoint configurations. This locally multiple Gō contact formalism is useful to investigate detailed pathways for conformational change. There can be ways to put a cooperativity in between these two limits: One can use the global multiple-Gō modeling for certain fragments of a protein [44,45]. Alternatively, one can design an arbitrary range of cooperativity in the conformational transition [46].

With the framework of constructing multiple-basins in hand, one can treat the ligand binding dynamic process in a few different manners. For a minimal modeling of conformational transition dynamics, one may simply change the relative free energy between two basins mimicking the ligand binding [47]. As a second level, one can introduce an implicit ligand binding interaction potential which effectively stabilizes the residues neighboring the ligand. Then, the ligand binding and unbinding is mimicked by truing on and off of this implicit ligand binding potential. This turning on and off process can be realized by Monte Carlo process, thus making the entire simulation as the hybrid-MC/MD simulation. Namely, the ligand binding and unbinding, or even chemical reactions, can be treated by MC, whereas protein conformational motions are treated by MD [48]. This scheme has been applied to modeling whole enzymatic dynamics of adenylate kinase as well as FO part of F-type ATP synthase [44,49]. As a third and more microscopic view, one can treat ligand as an explicit molecule represented by a few beads [42].

To mimic an entire cycle of conformational transition for huge molecular complexes, such as molecular machines, one may favor even a simpler approach; sudden switch of the reference (native) structures in Gō models. For example, a protein may have an open conformation A in the apo state and a closed conformation B in the ligand-bound state. For mimicking the ligand-binding induced conformational change, during a MD simulation with a Gō model based on the A structure, one can suddenly (or gradually) switch the potential surface to a new Gō model based on the B structure. After the sudden (or gradual) switch, the protein conformation smoothly relaxes from the A basin to the new structure B, which can be thought as a proxy to the conformational change. It should be noted that the sudden switch of potential places the protein at quite a high energy state in the new potential surface so that the initial part of the relaxation process tends to contain more artefacts. This “switching Gō model” was applied to a rotary molecular motor, F1-ATPase [50], molecular chaperone [51], and, more recently, for ATP-dependent chromatin remodeler [52].

Integrating with physical model

While the Gō model represents an ideal protein encoding the prefect funnel picture, real proteins do have some non-ideality. Functional restraints give rise to frustrations that are enriched near active sites as well as allosteric sites [30,53]. To model non-ideality on top of the funnel picture, it has been devised to integrate the Gō type bias with purely sequence-based and physico-chemical potentials. It should be noted that earlier off-lattice Cα Gō models such as [28] already took into account the sequence-based interaction to some extent. More thorough sequence-based interactions, i.e., non-native contact interactions, were also incorporated into Gō type modeling; for example, the AWSEM provides a spectrum of models from a purely sequence-based model for protein structure prediction to a Gō type model representing the perfect funnel landscape [54,55].

Within Cα models, one can approximate backbone hydrogen bonds and sidechain interactions using orientation-dependent potential functions. Hoang, T. X., et al. proposed to use three consecutive Cα positions to calculate orientations of backbone amide- and carbonyl- groups of the central residue so that the orientation-dependent backbone hydrogen bonds can be modeled as a function of 6 Cα coordinates, which was carefully calibrated later [56,57]. A similar approach was further developed to model the sidechain orientation in terms of three consecutive Cα positions [58,59]. This type of modeling is particularly useful to model intrinsically disorder proteins/regions that lack the native structure, as well as amyloid like higher-order structures. While these physico-chemical interactions alone cannot specific enough to fold a globular protein purely from sequence information, one can combine it with Gō type bias, whenever necessary.

Further spreading and future directions

While Gō models have originally been developed for and applied to proteins, one can extend the same idea to other macromolecules. The RNA folding and functional dynamics have been studied by some Gō models together with physico-chemical interactions [60-62]. For DNA, Gō type biases have been added to more physico-chemical interactions to stabilize B-type duplex form [63] although this bias was much weakened in more recent modeling. Notably, Gō type model has never been used for lipid, to my knowledge; lipid in biology stays in the liquid phase and thus it is not compatible with the idea of Gō model. In other words, the applicability of Gō models is linked to the availability of X-ray crystallographic structures; if a biomolecule crystalizes and its crystal diffracts X-ray well, that structure can be specific enough to rationalize the use of Gō models.

Gō models have also been used at the level of all-atom models [64]. Sanbonmatsu, Whitford, and Onuchic successfully applied it to reveal energy landscape of ribosome conformational dynamics and others [65]. The SMOG server provides a useful platform to run the all-atom, as well as coarse-grained, Gō models [66,67].

While Gō models were originally proposed as means to elucidate protein folding mechanisms, current and perhaps future usage of Gō models more often aims at fast conformation sampling. To characterize free energy landscape/pathways of some protein conformational change at atomic details, one often uses advanced sampling methods, such as Markov state models, and the string-method. These methods need initial and rough samples of relevant part of conformational space. Gō models are very efficient tools to prepare these samples.

Another interesting direction is structure-modeling; when high-resolution structures are available for most parts of a protein separately, but not the full-length of the protein, one can integrate them using a sort of Gō models [45]. Notably, most Gō models do not need information of Cartesian coordinates of the entire system. Instead, they need collections of many pairwise contact information. Thus, structure-information of many parts can be naturally merged into the full-length protein within Gō models

For many of these useful applications, when one uses Cα Gō model, one needs to reconstruct all-atom models that are compatible with the given Cα structure. This so-called back-mapping has been well developed and one can immediately reconstruct all-atom protein models given a reasonable Cα-model structure [68,69].

Publicly available tools for Gō models

There are several publicly available tools/servers to run Gō model simulations immediately although writing the code from scratch can also be straightforward. The SMOG server, a very useful web server (the SMOG can be used as a standalone program, as well), produces the necessary data for Gō model simulations with GROMACS [70], by which one can run MD [67,71]. It covers many of the variants described here, including Cα and all-atom Gō models, multiple-basin Gō model, as well as RNA Gō model. The MMTSB (Multiscale Modeling Tools for Structural Biology) tool set also offers a “Go Model Builder”, which provides necessary files for running Cα Gō model simulations with MD engines, such as CHARMM, GENESIS, and others [72-74]. The eSBMTools is a convenient python source-code package for various Gō model simulations, providing interface to run MD with GROMACS, somewhat similar to SMOG [75]. The CafeMol is a standalone software that implements a simple as well as a fine-tuned version of Cα Gō models, a multiple-basin Gō model, a RNA Gō model, together with an experiment-based coarse-grained DNA model and protein-DNA interactions [37]. The AWSEM-MD package offers a broad range of models from a purely physico-chemical and sequence-based interaction model, in one limit, to a Gō model in the other limit [55]. It uses a three-beads-per-amino-acid resolution. The AWSEM-MD produces an interface to run MD with the LAMMPS MD suite.