Efficient Traversal of Beta-Sheet Protein Folding Pathways Using Ensemble Models

Abstract

Molecular dynamics (MD) simulations can now predict ms-timescale folding processes of small proteins; however, this presently requires hundreds of thousands of CPU hours and is primarily applicable to short peptides with few long-range interactions. Larger and slower-folding proteins, such as many with extended β-sheet structure, would require orders of magnitude more time and computing resources. Furthermore, when the objective is to determine only which folding events are necessary and limiting, atomistic detail MD simulations can prove unnecessary. Here, we introduce the program tFolder as an efficient method for modelling the folding process of large β-sheet proteins using sequence data alone. To do so, we extend existing ensemble β-sheet prediction techniques, which permitted only a fixed anti-parallel β-barrel shape, with a method that predicts arbitrary β-strand/β-strand orientations and strand-order permutations. By accounting for all partial and final structural states, we can then model the transition from random coil to native state as a Markov process, using a master equation to simulate population dynamics of folding over time. Thus, all putative folding pathways can be energetically scored, including which transitions present the greatest barriers. Since correct folding pathway prediction is likely determined by the accuracy of contact prediction, we demonstrate the accuracy of tFolder to be comparable with state-of-the-art methods designed specifically for the contact prediction problem alone. We validate our method for dynamics prediction by applying it to the folding pathway of the well-studied Protein G. With relatively very little computation time, tFolder is able to reveal critical features of the folding pathways which were only previously observed through time-consuming MD simulations and experimental studies. Such a result greatly expands the number of proteins whose folding pathways can be studied, while the algorithmic integration of ensemble prediction with Markovian dynamics can be applied to many other problems.

1. Introduction

Protein folding and unfolding is a key mechanism used to control biological activity and molecule localization (Dobson, 2003). The simulation of folding pathways is thus helpful to decipher the cell behavior. Classical molecular dynamics (MD) methods (Karplus and McCammon, 2002) can produce reliable predictions but unfortunately the heavy computational load required by these techniques limits their application to inputs tens of amino acids long and prevents their application to large sequences (i.e., hundreds of amino acids). Recently, Faccioli et al. (2006) proposed an effective solution of the Fokker-Planck equation to compute dominant protein folding pathways, but the same size limitations remain.

The development of distributed computing technologies has dramatically extended the range of application of MD techniques. For instance, Pande and co-workers achieved a 1.5 millisecond folding simulation of a 39 residue protein NTL9 (Voelz et al., 2010). In spite of this achievement, this strategy still seems limited to small polypeptides (about 50 residues) and, more importantly, requires several months of parallel computing and typically thousands of GPU's.

In this article, we introduce a complete methodology to address these computational complexity limitations. Our approach aims to complement the range of techniques already offered. Unlike MD simulations that use an all-atom description of structures together with a fine-tuned energy force field, here we use a residue-level representation of the structure with a statistical residue contact potentials. This simplification enables us to sample intermediate structures and build a coarse-grained model of the energy landscape and subsequently simulate folding processes.

Since the seminal work of Levitt and Warshel (1975), it is widely acknowledged that simplified representations of protein structures and motions are required to circumvent computational limitations. A conceptual breakthrough came when Amato and co-workers applied motion planning techniques to the protein folding problem (Tapia et al., 2010; Amato and Song, 2002). The method is much faster than classical MD techniques and enables the study of the folding of large proteins. However, this approach does not predict structures; rather, it requires the three-dimensional (3D) structure of the native state to compute potential intermediate structures and unfolding pathways, on which the folding simulations are performed. It follows that the methodology cannot be applied to proteins with unknown structures and cannot be relied upon to study misfolding processes.

In fact, all the methods previously described face a difficulty common with MD: efficient sampling of the conformational landscape. MD algorithms explore the landscape through force-directed local search and progressive modification of the structure. However, the scalability and numerical efficiency when modeling large molecular structures remains problematic, limiting their application to small molecular systems. On the other hand, motion planning algorithms use a 3D structure of the native fold to predict distant structural intermediates. Accordingly, the accuracy of the method can suffer when intermediates sampled are far away from the native state. Recently, Hosur et al. (2011) have combined efficient motion planning techniques with machine-learning to model proteins as an ensemble, but this approach is effective only in the local neighborhood of the input structure.

Such obstacles have been addressed for RNA molecules by the development of structural ensemble prediction algorithms (McCaskill, 1990; Ding and Lawrence, 1999), and the derivation of a finely-tuned energy model based on experimental data (Turner and Mathews, 2010). Combined together, these techniques enable us to compute the RNA secondary structure energy landscapes and sample structures from sequence information alone. Wolfinger et al. (2004) further demonstrated how an RNA energy landscape can be constructed by connecting these samples together and estimating the transition rates between pairs of interconverting states. The resulting ordinary differential equation (ODE) system can be solved to predict and characterize RNA folding pathways. The method has since been improved to analyze the motion of large RNAs (Tang et al., 2008).

In this article, we propose to expand the methodology developed for RNAs to the more complicated case of proteins. First, we design an algorithm to sample the complete conformational landscape of large protein sequences given sequence data alone. Then, we use this sampling algorithm to build a coarse-grain representation of the energy landscape of a protein, from which we construct an ODE system modeling transition rates between folding intermediates that we solve to simulate protein folding.

We choose to address specifically β-sheet structures. The folding of these structures is particularly difficult to simulate. Indeed, β-sheets are stabilized by inter-strand residue interactions, and thus the folding and assembly of these structures is largely influenced by long-range interactions and global conformational rearrangements. For instance, Voeltz et al. (2010) recently showed that the rate-limiting step in the NTL9 fold was beta-sheet hairpin formation.

Since the original work of Mamitsuka and Abe (1994), several groups have proposed models to predict general β-sheets (Chiang et al., 2006; Kato et al., 2009; Tran et al., 2011). However, none of these methods are capable of computing ensembles of β-sheet structures (i.e., perform an exact enumeration of all β-structures without duplicates) and therefore cannot be used to sample the β-sheet energy landscape.

We recently introduced a structural ensemble predictor for transmembrane β-barrel (TMB) proteins (Waldispühl et al., 2008), continuing earlier work on molecular structure modeling (Waldispühl and Steyaert, 2005; Waldispühl et al., 2006). However, TMBs are a special case of β-sheets where each strand pairs with its two sequence neighbors via an anti-parallel interaction (except the “closing” pair which involves the first and last strands). Here, we expand these techniques to allow any β-strand organization in the β-sheet, with parallel and anti-parallel orientations, and enable the sampling of general β-sheets. This algorithm is implemented in the program tFolder.

We use tFolder to sample the β-sheet conformational landscape and build a coarse-grain model of the energy landscape. More specifically, we cluster protein configurations according to contact distance metrics, and associate each cluster with an intermediate folding state. We use the difference between the ensemble free energies of the clusters to compute the transition rates and build an ODE system that models the energy landscape. Finally, we solve this system to estimate the distribution of conformations over folding time.

This methodology reconciles the MD and motion planning approaches for studying folding pathways. Using tFolder, we are now able to simulate in a couple of minutes on a single desktop the folding of large proteins, and to predict the folding pathways (as well as possible misfolding pathways) of proteins with unknown structures. Thus, we are able to provide a broader range of applications, while offering computational efficiency comparable with motion planning techniques. Although we focus on β-sheet proteins, our method in principle could be extended to describe the folding pathways of a wider class of protein structures.

This article is organized as follows. In Section 2, we describe the tFolder algorithm and explain how we construct the coarse grained energy landscape model. Then, in Section 3, we benchmark our methods. First, we evaluate the accuracy of tFolder for simple inter-strand residue contact prediction and show that it performs comparably with more sophisticated techniques specifically designed for this task. Importantly, our contact predictions are not dependent on the separation between the residue indices, which means an improved “very” long-range contact prediction accuracy. Then, we illustrate the insights provided by our methods by analyzing the energy landscape of the extensively studied Protein G. We show that tFolder predicts the correct folding pathways, and interestingly, our simulation reveals a possible off-pathway structure. All these simulations can be performed on query sequences using our program tFolder, available at http://csb.cs.mcgill.ca/tFolder.

2. Methods

To predict realistic protein folding pathways, we exploit well-established ensemble prediction algorithms (Waldispühl et al., 2008) for their ability to accurately predict the energy scores of millions of feasible structural conformations from sequence alone. Our approach proceeds in two steps: (1) Given an arbitrary peptide sequence, we produce ensemble predictions of the energetic weight for all possible β-sheet structures and sub-structures, utilizing an enhancement to standard ensemble predictors which allows permutation. (2) Using each conformation's energetic score and metrics of conformational similarity, we derive the likelihood of dynamic state-to-state transitions and assemble a set of complete folding paths. In this way, we can identify and rank the most likely pathways from an unfolded conformation to a fully folded conformation based on predicted energy landscapes.

Modeling β-sheet ensembles

We model the set of all possible β-sheet conformations a peptide can attain using a statistical-mechanical framework. Conceptually, each structure is described by the set of residue/residue contacts that form hydrogen bonds between β-strand backbones, and is assigned a Boltzmann-distributed pseudo-energy, determined by the specific residues involved in contacts. To characterize the energetic landscape of this ensemble, a partition function Z can be calculated over all structural states \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$S = \{ 1 \ldots n \} $$\end{document} such that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}\begin{align*} Z = \mathop \sum_ { i = 1 } ^n e^ { - \frac { E_ { S_i } } { RT } } , \end{align*}\end{document}

with energies \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$E_{S_i}$$\end{document} , temperature T, and the Boltzmann constant R. For example, from this the relative abundance of a structure S_i can be easily derived: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}\begin{align*} p ( S_i ) = {\frac {e^ {- E_ { S_i }} } { Z } } . \end{align*}\end{document}

Our energy model is based on statistical potentials and follows directly from prior prediction tools that have been shown to be accurate (Waldispühl et al., 2008; Cowen et al., 2001). An energy E_i,j is given to each residue/residue pair within the β-sheet fold following E_i,j = −RT[log(p(i, j)) − Z_c], where Z_c is a statistical recentering constant and p(i,j) is the probability of these two residues appearing in a β-sheet environment, as observed across all non-sequence-homologous solved structures in the PDB (Waldispühl et al., 2008). Further, we assign separate probabilities based on the hydrophobicity of the environment on either face of a β-sheet.

A naive approach to computing the partition function would thus be to enumerate all possible structures and compute each structure's contribution to the sum individually. However, as was previously shown for the special case of anti-parallel β-strands in transmembrane β-barrel proteins, a much more efficient method exists using dynamic programming (Waldispühl et al., 2008). We have generalized this approach to enable the computation of arbitrary single β-sheet fold topologies.

Permutable β-templates

We introduce the concept of permutable β-templates to enable the calculation of the partition function of a β-sheet with arbitrary β-strand topologies. This extends existing ensemble prediction techniques by allowing any combination of parallel and anti-parallel β-strands to be including within a single β-sheet fold, and by removing any sequence dependency between β-strand /β-strand pairing partners. Prior methods supported only all-anti-parallel β-strands and required β-strand /β-strand interactions to be separated only by coil (and not other strands) (Waldispühl et al., 2008).

To efficiently encode these generic shapes, each strand is labeled \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\{ 1 \ldots n \} $$\end{document} to allow a stepwise permutation through β-strand ordering, and a signed permutation is defined such that each β-strand is assigned to be parallel or anti-parallel relative to the first strand in the sheet (Fig. 1). Algorithmically, tFolder is capable of constructing a dynamic program over all such permutations to calculate the partition function. In practice, since such an encoding can result in unrealistic combinations of β-strand/β-strand pairings (such as if β-strands 1 and 4 had too short a coil between them in Fig. 1), we impose that valid foldings must satisfy steric and biologically derived constraints. These include a minimum and maximum β-strand length, maximum shear between neighboring β-strands (the amount of inclination that causes the β-sheet to deviate from a perfect rectangle), and minimum inter-strand loop size. These constraints serve to limit the exploration of unrealistic conformations, minimizing excess computation and allowing directed investigation into specific motifs.

FIG. 1.

An illustration of how a permutable β-template can be encoded as a signed permutation. The permutation lists the strands in the order that they occur in the sheet, with the sign indicating whether the strand is parallel (+) or anti-parallel (−) to the first strand.

The energy of a structure with n strands, can be recursively defined as E(S_n) = E(S_n₋₁) + Pairing(s_n₋₁,s_n), where E(S_n₋₁) is the interaction energy between the first n − 1 strands, and Pairing(s_n₋₁,s_n) is the energy of the pairing of strand n − 1 with strand n (Fig. 2a). tFolder exploits the shared structure between instances in the ensemble by computing this recursion using a dynamic programming algorithm. The result of each recursive call is stored in a table indexed by the parameters of the call. Subsequent recursive calls made with the same parameters perform a table lookup instead of re-computing the value of the recursion.

FIG. 2.

(a) Illustrates how the energy function of a β-sheet can be recursively defined as the sum of the contribution of the last two strands with the contribution of the remaining structure. (b) Indicates the indices used to store the energies of intermediate structures for the recursion.

For a sheet of n strands, the table has n rows, where the kth row has entries corresponding to valid configurations of the first k strands. For the kth strand, these configurations are partitioned by the location of four indices k₁, k₂, k₃, k₄, which denote the boundaries of the region occupied by the k strands (Fig. 2b). To begin, the algorithm enumerates all possible positions of the first two strands, and for each stores the strand pair interaction energy in entry E_2₁_2₂_2₃_2₄ of the table. For each subsequent strand k, the value of E_k_1k2k3k4 is computed as: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}\begin{align*} E_{k_{1}k_{2}k_{3}k_{4}} = \mathop \sum_{i_{1}i_{2}i_{3}i_{4}} E_{i_{1}i_{2}i_{3}i_{4}} + Pairing ( i , k ) , \end{align*}\end{document}

where i₁, i₂, i₃, i₄ are enumerated for all valid settings for the boundaries of the preceding strands, given the boundaries of the kth strand. Once the recursion has filled the table, the partition function Z is calculated by summing over all possible settings of n₁, n₂, n₃, n₄: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}\begin{align*} Z = \exp \left( - \mathop \sum_{n_{1}n_{2}n_{3}n_{4}} E_{n_{1}n_{2}n_{3}n_{4}} \right). \end{align*}\end{document}

The table constructed to calculate the partition function can be used to sample the distribution of configurations of a given topology, utilizing the approach established by Ding and Lawrence (2003) for RNA secondary structure, and successfully applied previously by Waldispühl et al. (2006) to sample conformations of β-barrel proteins. To do this, we perform a traceback through the table and, at ith step, sample the indices within which the first i strands are contained, according to the Boltzmann representation of these i-stranded structures (Fig. 3).

FIG. 3.

Illustration of of how the sampling procedure performs a traceback through the table, over the indices of intermediate structures. During each step of the sampling procedure, the location of a single strand is sampled from the region indicated by the vertical bars. The triangles denote the location of the strand sampled during the previous step.

2.1. Predicting folding dynamics

Conceptually, we model the folding process as a path through a graph of varyingly folded conformations of a protein. In this graph, different protein conformations are represented as states, and two states that inter-convert in a folding pathway are connected by an edge, analogous to work with RNA described previously (Wolfinger et al., 2004). The tFolder algorithm provides a means to efficiently sample the energetically accessible conformations that make up the states of this graph. We further propose a means to determine the connectivity between states and demonstrate how this can be applied to calculate the dynamics of the folding process.

Since we do not know the final structure, we begin by sampling configurations from all possible permutations of β-sheet topology, as described above. For every pair of states, we add an edge between two states if (1) the states have compatible topologies, and further, (2) the states show structural similarity.

Two templates are compatible if they are identical to each other, modulo the addition or removal of a single strand pairing. This operation can result in the growth of a core structure, or the nucleation of an independent strand pair (Fig. 4). Note that the requirements for satisfying the second criterion of structural similarity depends on the metric used to estimate structural similarity between two conformations. In practice, we use a contact based metric and deem two structures to be structurally similar if the metric is below the transition threshold.

FIG. 4.

The topologies that are compatible with a given state (shaded gray) result from the addition of a single pairing between strands (dashed box). The “+” indicates that there is no pairing between the gray structure and the white strand pair.

Given the graph constructed according to these two criteria, the change in the probability of the system being in state i at time t is calculated from the total flux into and out of state i, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}\begin{align*} \frac { dp_ { i } } { dt } = \mathop \sum _ { j \epsilon X } r_ { ij } p_ { j } ( t ) , \end{align*}\end{document}

where p_i is the probability of state i, X is the state space, and r_ij is the rate of transition from state i to state j. Given that two states are connected in the graph, the rate at which two states inter-convert is proportional to the difference between free energies of the states (ΔG); the system tends toward energetically favorable states. We calculate the transition rate r_ij between states i and j using the Kawasaki rule (with parameter r₀ to scale the time dimension): \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}\begin{align*} r_{ij} = r_{0} \exp \left( - \Delta G_{ij} / 2 RT \right). \end{align*}\end{document}

The dynamics of the system are calculated by treating the folding process as a continuous time discrete state Markov process. Given the matrix of folding rates R, where R_ij = r_ij and initial state density p (0), the distribution over states p (t) of the system at time t is given by the explicit solution to the system of linear differential equations, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}\begin{align*} \textbf{\textit{p}} ( t ) = \exp ( Rt ) \textbf{\textit{p}} ( 0 ). \end{align*}\end{document}

Since we sample hundreds of states from each β-strand topology, we partition the state space into macro states using clustering, in order to work with a tractably sized system. Under this approximation, we consider two clusters the graph to be connected if the minimum distance between any two states from each cluster is below the transition threshold. We define the ensemble free energy difference ΔG_ij between two macrostates i and j by summing over the states from which they are composed. \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}\begin{align*} \Delta G_{ij} = E ( \chi_{i} ) - E ( \chi_{j} ) = \mathop \sum_{x \epsilon \chi_{i}} E ( x ) - \mathop \sum_{x \epsilon \chi_{j}}E ( x ). \end{align*}\end{document}

Although this approximation lessens the computational burden, it represents a trade off. The granularity achievable by our simulation is at the level of the macrostates. Note, energy barriers are not explicitly incorporated into the model, since entire β-strands are either added or removed between states without partially-formed intermediates.

3. Results

3.1. Evaluation of contact prediction

To evaluate the contact prediction performance of tFolder, we tested it using a 16-protein benchmark.¹ Proteins were selected from the Protein Data Bank which had 4–6 β-strands with dominantly beta structure, 100 residues or smaller, and sequence identity less than 30%. To ensure a fair comparison with the BETApro and SVMcon, we removed from this set those proteins that were used in the training sets of these methods. From each of these 16 proteins, the β-topology was extracted and used as input for tFolder, along with the amino-acid sequence and fixed strand length of 5–8 residues. Since predicting folding dynamics involves a permutation over all β-topologies, this demonstrates the expected accuracy of each folding state along the pathway. We sample 500 configurations of each protein, and use these ensembles to compute a stochastic contact map and distribution of strand locations (Fig. 5a, and b). The contact map represents the probability of observing a given contact, and predicted contacts are the set of all contacts with probability above a threshold value t. The selection of t influences the measured performance (Fig. 5c), so to objectively set the threshold, we chose a t that maximizes the F-measure. We evaluated the quality of our contact maps based on the Accuracy \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\left( \frac { no. \ of \ correctly \ predicted \ contacts } { no. \ of \ predicted \ contacts } \right)$$\end{document} , Coverage \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\left( \frac { no. \ of \ correctly \ predicted \ contacts } { no. \ of \ observed \ contacts \ contacts } \right)$$\end{document} , and F-measure \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\left( \frac { 2 \cdot Accuracy \cdot Coverage } { Accuracy + Coverage } \right)$$\end{document} of our predictions. We calculated these measures in terms of β-contacts, which we defined as residues located within β-strands less than 8Å apart (between C_α atoms) in the PDB structure. A summary of tFolder performance on Protein G, as well as average performance on the 16 protein dataset, is presented in Table 1. Here we distinguish between results for long range contacts, greater than 0, 12, or 24 residues apart. Thus, tFolder maintains reasonable predictive accuracy even with large contact separations.

FIG. 5.

Summary of the distribution of structures predicted by tFolder for Protein G. (a) Contact probability predicted by tFolder between all pairs of amino acids. Green squares indicate contacts predicted by tFolder, whereas red squares represent pairs of amino acids that are less than 8Å from each other in the observed structure. A higher intensity of green indicates a higher predicted probability of the contact, and yellow squares are an indication of agreement between prediction and observed contacts. (b) Probability of the location of each strand, computed from the ensemble of sampled structures. The bars at the top of the plot indicate the location of the strands from the experimentally determined structure. (c) Relationship between the threshold used to determine a contact from the contact map, and the values of the three metrics Accuracy, Coverage, and F-measure. The threshold can be set to maximize the value of the F-measure, representing a reasonable trade-off between Coverage and Accuracy.

Table 1.

Performance of tFolder for Contact Prediction Is Evaluated Based on the Accuracy, Coverage, and F-Measure of Experimentally Observed Contacts

	Protein G			16 protein benchmark
	≥0	≥12	≥24	≥0	≥12	≥24
F-measure	0.36	0.37	0.45	0.25	0.27	0.23
Accuracy	0.34	0.33	0.41	0.25	0.27	0.28
Coverage	0.39	0.42	0.50	0.28	0.32	0.25

These performance metrics are reported for contacts that are more than 0, 12, and 24 residues apart, showing that tFolder maintains reasonable predictive accuracy even with large contact separations.

In order to evaluate the performance of tFolder with respect to other approaches for contact prediction, we present in Table 2 a comparison with two leading contact prediction algorithms, SVMcon and BETApro. It can be seen that tFolder is able to perform comparably at the task of contact prediction, and performs particularly well for contacts with sequence separation greater than 24 residues. Although these methods sometimes outperform tFolder, it is important to note that the accuracy of tFolder is less sensitive to the distance of contact separation. Since critical protein folding steps can involve both short-range and long-range β-sheet contacts, it is especially important for long-range contacts to be predicted correctly to allow an accurate folding pathways to be reconstructed.

Table 2.

Comparison of the Performance of tFolder Contact Prediction with Contact Prediction Algorithms SVMcon and BETApro

	≥12			≥24
Method	F-measure	Accuracy	Coverage	F-measure	Accuracy	Coverage
tFolder	0.27	0.27	0.32	0.23	0.28	0.25
BETApro	0.22	0.40	0.16	0.05	0.14	0.07
SVMcon	0.32	0.31	0.50	0.24	0.21	0.44

The methods are evaluated based on their ability to perform contact prediction for contacts greater than 12 and 24 residue separation respectively.

3.2. Predicting the folding pathways of the B1 domain of Protein G

To demonstrate the efficacy of our techniques for predicting protein folding pathways, we reconstruct the folding landscape of the B1 domain of Protein G—a well-studied protein for which the pathway has been elucidated through many experimental studies and MD simulations. To do this, all possible permutations of a 4-strand β-sheet topology were sampled and clustered. For each of these sets of structures, the cluster with the highest probability of being observed was selected to be representative of each topology.

The graph of the folding pathway was constructed by considering all pairs of clusters. If the minimum distance beween two clusters was less than a predetermined transition threshold, we considered that there was exchange between the two states. We tried several metrics, including segment overlap, mountain metric, and a contact based metric (Zemla et al., 1999; Moulton et al., 2000), selecting the contact-based metric, because it performed best empirically. The resulting graph of protein conformations is illustrated in Figure 6a. Inspection of this graph, along with the folding dynamics computed from this graph in Figure 6b, reveals folding intermediates consistent with those previously reported by Song et al. (2003). It should also be noted that although we compute other configurations of the sequence that are energetically favorable (faded states in Fig. 6a), they are not predicted to form because they are unreachable from the unfolded state. Interestingly, a four-stranded off-pathway structure is predicted to form, which has not been observed previously. Furthermore, our results agree with the work of Hubner et al., (2004), who show that the anti-parallel beta-hairpin, predicted to form an interaction between residues 39–44 and 50–55, center around known nucleation points W43, Y50, and F54.

FIG. 6.

(a) Graph of the folding landscape of Protein G predicted by tFolder is illustrated above. The gray shaded region indicates the states predicted to be reachable from the unfolded state. The dark arrows indicate transitions between states, and the size of the arrow indicates the favored direction of transition along each edge. Faded arrows are drawn between states that have compatible topologies but do not reach the transition threshold. The size of each state indicates its relative representation at equilibrium. The faded structures indicate states that are unreachable from the unfolded state. (b) Folding dynamics of Protein G shows how the probability of observing any of the reachable states changes over the time the protein folds. Each line is annotated with an image of the state it represents.

3.3. Algorithm running time

The computational bottleneck of our approach is the computation of the partition function of a template. The primary factors influencing this calculation are the length of sequence and the number of strands in the β-topology (the depth of the recursion). The partition function for sequences between 40–130 residues and 4–6 strands was calculated using a single 2.66-GHz processor with 512 MB of RAM. The effect of these two parameters on the computation time is depicted in Figure 7 below. Further, computing the parition function across multiple β-templates is trivially parallelizable. The ability to formulate quick, coarse-grained predictions in a matter of minutes, rather than days of atomistic-detail simulation, is a fundamental benefit of our technique.

FIG. 7.

The time required to compute the partition function increases with increasing size of amino acid sequence, and number of strands. The time was computed by averaging over n = 3 trials, for sequences of 40–130 residues in length, with 4–6 strands.

4. Discussion

We present tFolder, a novel approach for quickly predicting protein folding pathways through the accurate prediction of the conformational landscape of arbitrary β-sheet proteins. What distinguishes tFolder from other computational approaches that attempt to probe protein folding processes is that tFolder does not require vast computational resources; in fact, it can be ran on a single personal computer. To achieve this performance we use a simplified model for protein folding, allowing us to very rapidly compute a coarse-grained picture of the folding of a protein from sequence information alone. This contrasts with methods that attempt to determine folding mechanisms by trying to unfold proteins from their native state. Such methods require the a priori knowledge of the native structure, and as such are not applicable to study protein sequences with unknown structures. When computing protein folding pathways, our method explores all possible β-sheet configurations, and thus does not face this limitation. Interestingly, this independence from known structures could provide insights into off-pathway kinetics, such as the aggregation of proteins into amyloid structures.

Although tFolder only predicts coarse folding pathway transitions in β-sheet proteins, its strength lies in its ability to quickly separate conformational transitions that are critical to folding from those transitions that could simply result from minor structural fluctuations. This complements the use of MD simulations as the MD can be used to explore the nuanced structural interactions that certainly occur near a transition highlighted by tFolder. Further, although we are able to produce good results using a fairly simplistic energy model, a more complicated formulation, such as one including entropic forces, would clearly improve tFolder's analysis. More advanced heuristics also exist (Tang et al., 2008), which more efficiently extract folding pathway information, which could be applied to tFolder.

Understanding the folding dynamics of β-sheet proteins, especially which β-strand contacts drive folding and conformational stability, could help create better models of hierarchical folding, protein aggregation, and evolutionary pressure. Significant overlap likely exists between many proteins' folding pathways to even permit a classification of common transition elements (Fulton et al., 2005); however, creating such a database would only be possible with sufficiently fast and accurate algorithms. tFolder takes a step toward this end by demonstrating techniques for efficiently predicting ensembles of arbitrary β-sheet proteins, and for combining these predictions to construct accurate protien folding transition landscapes.

5. Appendix

Computing the partition function

For a β-sheet with k strands, and residues labeled \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$1 \ldots n$$\end{document} , its topology can be represented as a vector p containing a permutation of the numbers \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\{ 1 \ldots k \} $$\end{document} . The partition function of this sheet is computed by a recursion of depth k, where at each level of the recursion the location of the strands p_i is summed over. In order to calculate this pairing energy, one must maintain the location of the last strand added to the structure, as well as the locations in the sequence that are already occupied by strands (to ensure that no strands are overlapping). The region already occupied by strands is represented by a vector of four indices {i₁, i₂, j₁, j₂}, indicating the boundaries of this region. This scheme allows one to track the boundaries of four boundaries, meaning it cannot be applied to all possible permutations. For instance if p = {1, 3, 5, 2, 4}, one would need to maintain the location of six boundary locations at the third level of the recursion. Nevertheless, keeping four boundaries seems sufficient in practice for the majority of observed sheet topologies, and if desired a more general structure could be developed at the price of carrying more information between levels of the recursion.

Depending on the permutation, the location of ith strand will occur in one of three regions relative to the boundaries of the occupied structure:

1. In the region of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$1 \ldots i_1$$\end{document}

2. In the region of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$j_1 \ldots i_2$$\end{document}

3. In the region of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$j_2 \ldots n$$\end{document}

Once the energy of the pairing has been computed, the boundaries of the occupied region are updated depending on which region the strand was added to, and whether further strands will be added in the gaps between strands. For example, in the first column, if there remains a strand to be added in the region \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$j_1 \ldots i_2$$\end{document} according to p, the boundary will be updated to {i₃, j₁, i₂, j₂}, otherwise the bounds will be {i₃, j₃, i₁, j₂}. In this manner the energy of a sheet can be recursively decomposed into a sum of pairing energies between strands. Thus, by selecting the proper sequence of decompositions, and enumerating the possible strand locations at each step, one is able calculate the partition function of the desired β-sheet. For clarity, illustrations of these decompositions and accompanying equations describing the energy of a partial structure are illustrated in Figures 8 and 9.

FIG. 8.

Illustration of the decompositions applied by tFolder used to compute the partition function. The gray arrow represents the addition of the ith strand, and the white arrow is the location of the previous strand, which must be maintained to compute the pairing energy. The solid arcs between (i₁, j₁) and (i₂, j₂) represent the boundaries of the previously occupied structure. The dashed arc filled arc represent the updated boundaries after adding the ith strand.

FIG. 9.

Equations used in the recursive calculation of the folding energy and partition function. The three groups separated by bars represent the three columns in Figure 8.

Footnotes

Disclosure Statement

No competing financial interests exist.

1

Complete benchmark results available at the tFolder website.

References

Amato

N.M.

, Song

2002. Using motion planning to study protein folding pathways. J. Comput. Biol., 9:149–168.

Chiang

, Joshi

A.K.

, Searls

D.B.

2006. Grammatical representations of macromolecular structure. J. Comput. Biol., 13:1077–1100.

Cowen

, Bradley

, Menke

et al. 2001. Predicting the beta-helix fold from protein sequence data. J. Comput. Biol., 8:261–276.

Ding

, Lawrence

C.E.

1999. A bayesian statistical algorithm for RNA secondary structure prediction. Comput. Chem., 23:387–400.

Ding

, Lawrence

C.E.

2003. A statistical sampling algorithm for RNA secondary structure prediction. Nucleic Acids Res., 31:7280–7301.

Dobson

C.M.

2003. Protein folding and misfolding. Nature, 426:884–890.

Faccioli

, Sega

, Pederiva

et al. 2006. Dominant pathways in protein folding. Phys. Rev. Lett., 97:108101.

Fulton

K.F.

, Devlin

G.L.

, Jodun

R.A.

et al. 2005. PFD: a database for the investigation of protein folding kinetics and stability. Nucleic Acids Res., 33:D279–D283.

Hosur

, Singh

, Berger

2011. Sparse estimation for structural variability. Algorithms Mol. Biol.(in press).

10.

Hubner

I.A.

, Shimada

, Shakhnovich

E.I.

2004. Commitment and nucleation in the protein G transition state. J. Mol. Biol., 336:745–761.

11.

Karplus

, McCammon

J.A.

2002. Molecular dynamics simulations of biomolecules. Nat. Struct. Biol., 9:646–652.

12.

Kato

, Akutsu

, Seki

2009. Dynamic programming algorithms and grammatical modeling for protein beta-sheet prediction. J. Comput. Biol., 16:945–957.

13.

Levitt

, Warshel

1975. Computer simulation of protein folding. Nature, 253:694–698.

14.

Mamitsuka

, Abe

1994. Predicting location and structure of beta-sheet regions using stochastic tree grammars. Proc. ISMB, 276–284.

15.

McCaskill

1990. The equilibrium partition function and base pair binding probabilities for RNA secondary structure. Biopolymers, 29:1105–1119.

16.

Moulton

, Zuker

, Steel

et al. 2000. Metrics on RNA secondary structures. J. Comput. Biol., 7:277–292.

17.

Song

, Thomas

, Dill

K.A.

et al. 2003. A path planning-based study of protein folding with a case study of hairpin formation in protein G and L. Proc. Pac. Symp. Biocomput., 240–251.

18.

Tang

, Thomas

, Tapia

et al. 2008. Simulating RNA folding kinetics on approximated energy landscapes. J. Mol. Biol., 381:1055–1067.

19.

Tapia

, Thomas

, Amato

N.M.

2010. A motion planning approach to studying molecular motions. Commun. Inform. Syst., 10:53–68.

20.

Tran

V.D.

, Chassignet

, Sheikh

et al. 2011. Energy-based classification and structure prediction of transmembrane beta-barrel proteins. Proc. First IEEE ICCABS.

21.

Turner

D.H.

, Mathews

D.H.

2010. NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure. Nucleic Acids Res., 38:D280–D282.

22.

Voelz

V.A.

, Bowman

G.R.

, Beauchamp

et al. 2010. Molecular simulation of ab initio protein folding for a millisecond folder NTL9(1–39) J. Am. Chem. Soc., 132:1526–1528.

23.

Waldispühl

, Berger

, Clote

et al. 2006. Predicting transmembrane beta-barrels and interstrand residue interactions from sequence. Proteins, 65:61–74.

24.

Waldispühl

, O'Donnell

C.W.

, Devadas

et al. 2008. Modeling ensembles of transmembrane beta-barrel proteins. Proteins, 71:1097–1112.

25.

Waldispühl

, Steyaert

J.-M.

2005. Modeling and predicting all-alpha transmembrane proteins including helix-helix pairing. Theor. Comput. Sci., 335:67–92.

26.

Wolfinger

M.T.

, Andreas Svrcek-Seiler

W.A.

, Flamm

et al. 2004. Efficient computation of RNA folding dynamics. J. Phys. A Math. Gen., 37.

27.

Zemla

, Venclovas

, Fidelis

et al. 1999. A modified definition of SOV, a segment-based measure for protein secondary structure prediction assessment. Proteins, 34:220–223.