Distinguishing biomolecular pathways and metastable states
Protein folding occurs in a high dimensional phase space, and the representation of the associated energy landscape is non-trivial. A widely applied approach to studying folding landscapes is to describe the dynamics along a small number of reaction coordinates. However, other strategies involve more elaborate analysis of the complex phase space. There have been many attempts to obtain a more detailed representation of all available conformations for a given system. In this work, we address this problem using a metric based on internal distances between amino acids to describe the differences between any two conformations. Using an effective projection method, we are able to go beyond the usual one-dimensional representation and thus visualize landscapes in two dimensions. We refer to this method as Energy Landscape Visualization Method (ELViM). We have applied this methodology using C α structure-based model (SBM) 1 to study the folding of two well-known proteins: SH3 domain and Protein-A. Our visualization method yields a detailed description of the folding process, making possible the identification of transition state regions, and establishing the paths that lead to the native state. For SH3, we have analyzed structural differences in the distribution of folding routes. The competition between the native and mirror structures in protein A is discussed. Finally the method is applied to discuss the conformational changes of the Elongation Factor Thermally unstable (EF-Tu) and its transition states. Distinct features of the ELViM are that it does not require or assume an existence of a reaction coordinate, and it does not require kinetic information regarding the investigated system.
Introduction<!>Results<!>Detecting Parallel Routes During Protein Folding<!>Folding in a Multi-basin Landscape<!>Capturing competing pathways during domain rearrangements<!>Discussion<!>Methods<!>I. Generating a data ensemble<!>II. Dissimilarity matrix<!>III. Preliminary data clustering<!>IV. Multidimensional reduction
<p>Protein folding occurs in a high-dimensional phase space, and it involves a conformational search across a broad range of structurally-distinct states. Within this vast space, the energy landscape is shaped such that a protein can reliably fold and perform its various functions. [1][2][3] Through the course of evolution, proteins have achieved this ability by selecting for energy landscapes that are dominated by the contacts formed in the native conformation. Thus, a long-term goal of the field has been to understand the character of folding dynamics, in the context of "minimally-frustrated", or "funnel-like" energy landscapes. 4,5 However, even though with this general character of the landscape, the precise dynamics results from an interplay of energetic, entropic and topological aspects. Thus investigating the balance between these factors has inspired the development of a large number of energetically "smooth" models for studying protein dynamics, 6 which includes a broad range of native-centric "structure-based" models. While these models are simplified, comparison of simulations 7,8 and experimental measurements 9,10 has shown that these models are able to capture many aspect of the folding process, including folding rates, Φ-values, landscape roughness, and topological impacts on folding mechanisms. [11][12][13] It is important to recognize that energy landscape techniques are not restricted to the study of protein folding. Rather, they may also be used to study functional dynamics, which can involve large-scale conformational changes, 14 motor-like energy transfer 15 and assembly. 16 In addition, this strategy can be used to investigate a variety of other biomolecular systems, including RNAs, 17 DNAs, 18 and complex assemblies, such as the ribosome 19 and chromatin. 20 Accordingly, these models represent a means for characterizing the dynamics of any number of large-scale and complex biomolecular processes.</p><p>Even though folding and functional dynamics are intrinsically multi-dimensional, it is often possible to accurately describe the kinetic and thermodynamic properties in terms of a few key quantities. For example, structural parameters, such as radius of gyration and fraction of native contacts (Q), are commonly used to measure protein folding dynamics. 21 While these low-dimensional representations can capture the dominant kinetic properties of a system, 22,23 and topological features of transition state ensembles, 11 applying putative coordinates has the potential to mask the richness of the dynamics. However, techniques are available for systematically identifying suitable reaction coordinates. For example, timecorrelation analysis 24 can be used to define classes of collective variable and then identify which describe the slowest motions. Transition path analysis may also be used to determine which coordinates capture an underlying free-energy barrier, 25 though exhaustive coordinate comparisons can become prohibitively expensive. One limitation with these techniques is that these methods require general definitions of coordinates be provided a priori. Accordingly, there is still a need for automated techniques that are capable of rapidly mapping the broad character of an arbitrary landscape.</p><p>An alternate strategy for mapping the landscape is to construct connectivity maps that are based on energetics or kinetics. For example, one may define Markov State Models (MSMs), where the system is described in terms of large numbers of relatively long-lived states. [26][27][28] Another approach is to identify local minima individually in order to graph the high-dimensional space. 29,30 For both of these techniques, more accurate representations are obtained when studying rough landscapes that possess many well-defined free-energy minima. However, when applying minimally-frustrated models, the minima tend to be characterized by broad and shallow minima, making the utility of these approaches unclear.</p><p>In the current study, we describe a method for visualizing energy landscapes, called the Energy Landscape Visualization Method (ELViM). While this technique has previously been applied to lattice models of proteins, 31 the current work extends this technique to study off-lattice models for folding and conformational rearrangements ELViM seeks to convey a meaningful 2D (and/or 3D) projection of the landscape that is based on an internal metric to assign distances between sampled conformations. In contrast to other methods, ELViM does not require intuitive designation of coordinate classes, nor does it require the sampling of long-lived minima, making it technique for studying smoother landscapes. To illustrate the applicability of this approach, we apply the method to three classes of conformational processes in proteins. First, we describe the folding of a well-studied protein: SH3. Next, we apply the method to study the folding of protein A, which can fold to two competing configurations: "native" and "mirror". Finally, we apply ELViM to describe large-scale (∼60 Å) domain rearrangements in elongation factor Tu (EF-Tu). Together, these applications demonstrate how ELViM enables broadly-applicable analysis that can distinguish the rich characteristics of biomolecular dynamics.</p><!><p>One goal of the present study is to develop and apply a method that is capable of identifying parallel pathways in nominally-smooth landscapes. As described in the introduction, simplified structure-based models can provide physically-meaningful descriptions of dynamics that are often consistent with experimental results. However, the challenge of identifying paths that are associated with broadly-defined ensembles is a difficult task. To approach this, we apply ELVIM, an analysis approach inspired by our earlier studies of lattice models for proteins. As described below, we find that this technique allows one to infer details of energy landscapes in a variety of contexts, including protein folding, folding to competing structures and large-scale rearrangements in a multi-domain protein.</p><p>An essential aspect of this approach is to have an appropriate metric to represent effective distances between pairs of conformations in a sampled data set. In this study, such metric corresponds to the internal distance introduced by the Wolynes group. 32 Using this distance measure, a multidimensional reduction is performed to obtain a 2D landscape representation, which provides easily interpretable descriptions of complex systems. In the present study, we restricted the metric to consider only structural information of each conformation, though subsequent extensions could include kinetic weighting. As we show below, even though kinetic information is not considered, the method still allows one to capture freeenergy barriers, transition states and parallel paths associated with folding and functional rearrangements.</p><!><p>When describing folding, the "mechanism" consists of a distribution of parallel transition paths. 33 There are cases in which visual inspection can suggests appropriate effective coordinates for identifying these parallel folding processes. However, it is desirable to establish protocols and methods that can identify parallel routes without prior knowledge of the system. One example of parallel folding paths can be found in the context of the SH3 domain of tyrosine kinase c-Src. The folding mechanism of SH3 has been well characterized, both experimentally and computationally, [34][35][36][37] making it an ideal benchmark for evaluating the robustness of the ELViM methodology. SH3 is a 57-residue protein that forms a five stranded beta barrel (PDB ID: 1FMK, 38 Fig. 1b). To study the folding of SH3, we apply molecular dynamics simulations using a C-α structure-based (Gō-like) model. 11 When using this model, the folding dynamics are two-state in character, and there are sharp transitions between the folded (high Q) and unfolded (low Q) ensembles (Fig. 1a). Similar dynamics is also observed when using the RMSD from the folded conformation (Fig. 1a). Further, when describing the dynamics in terms of both Q and RMSD (Fig. 1b), there does not appear to be an obvious When we use ELViM to project SH3 conformations into a 2D map, a rich picture of the underlying landscape is revealed, where the unfolded and folded ensembles are connected through a number of transition regions (Fig. 2a). One may also represent projection as a 2D free energy, where ∆F = − log(n i /N ), n i is a local density of points, and N is the total number of conformations. This reveals three major paths (A, B and C) with similar transition barrier heights (∼3 k B T) separating the unfolded and folded ensembles (Fig. 2b).</p><p>Previous efforts by Lammert and co-workers 36 showed that the so-called 'early' ('late') group of contacts are preferably formed before (after) the transition state is reached (Fig. 3d). To describe the dominant folding path, they introduced a reaction coordinate Q path , which accounts for the difference between the number of 'early' and 'late' contacts that are formed. The main transition paths shown in Fig. 2b (paths A,B, and C) are all in agreement with the initial formation of early contacts (Fig. 3d). Folding through the reverse mode (i.e.</p><p>forming "late" contacts first) is associated with a much larger free-energy barrier, which indicates that such a process is not likely to be kinetically competitive (see Fig. 3 and S2).</p><p>While all three transition state regions (A,B, and C) are associated with formation of early group contacts, ELViM reveals distinct differences in these parallel paths. For example, as can be seen by a representative structure of each region, shown in the lower part of each contact map (Fig. 3.a, b, and c). These contact maps can give a weighted transition path distribution, which is more meaningful than the average transition-state ensemble (TSE) representation using Q, which eventually is not even the best reaction coordinate.</p><p>ELViM takes into account only geometric distances between conformations. If one considers the potential energy associated with each structure, a third dimension is included in the landscape description. The SH3 3D representation presents indeed a funnel like structure, in which the transition state conformations are placed asymmetrically around the bottom of the funnel. This representation can help to visualize the favorable folding pathways, connecting them to the different transition states, see Fig. S3, and Video S1.</p><!><p>Mirror-like protein folds can form when there are energetically competitive structures accessible that are compatible with the native contact network. 39,40 As an example, using molecular dynamics simulations, Noel et al. 41 concluded that the native fold and its mirror image have similar stabilities, making them thermodynamically competitive.</p><p>To investigate how the interplay of competitive conformations may be described by 2D ELViM projections, we performed simulations using a dual-basin structured based model (model described in 41 ), where the native and mirror-like conformations represent free-energy minima.</p><p>To determine whether ELViM can distinguish between the competing folding routes, the folding landscape of Protein A (Fig. 4) is colored according to</p><p>Native or Mirror) is defined as the sum of all tertiary contacts present in each contact map (coord defined in Ref. 41 ). The region in blue represents conformations similar to mirrored structures, while the red region describes native-like structures. Given the parameters of the dual basin model, 41 simulation was conducted at temperature below T f , in which sampling of the Native and Mirror conformations are better probed, though completely unfolded conformations also occur (see Methods and SI). Fig. 4a shows that native-like conformations occur more frequently than the mirror-like ones. Unfolded conformations are all on the outer rim of the figure (see Fig. S4) A contour plot of F (Fig. 4b) reveals three clear minima: Native (N), a native-like intermediate (I) and the mirror image (M). The emergence of a native-like intermediate is unexpected, and it was not reported in earlier studies that used conventional low-dimensional reaction coordinates (Q, the Q N − Q M , and RMSDs). 41 This suggest a new mechanism in which the mirror and native states inter-convert into each other.</p><p>When one inspects representative structures of the I state, one can envision a process by which the Helix 3 (blue Helix) is displaced relative to a nominally-folded Helix 1/2 bundle (red/white bundle) (Fig. 4). When we take into account the potential energies of each state, the 3D visualization of Protein A energy landscape gives insight also on the relative energy differences between mirror and native configurations (Fig. S5 and Video S2). In this example we have seen how this reaction-coordinate free landscape analysis helped to identify an intermediate state, which may help to design probes to detect them experimentally.</p><!><p>Elongation factor-Tu (EF-Tu) is a three-domain protein that is responsible for delivery of aminoacyl-tRNA molecules to the ribosome. [42][43][44] In order to accomplish this task, EF-Tu binds aa-tRNA and GTP to form the ternary complex. During initial association with the ribosome, the aa-tRNA forms base pairs with the ribosome-bound mRNA codon. EF-Tu then hydrolyzes GTP and undergoes a large-scale domain rearrangement process, resulting in the release of the aa-tRNA molecule. Overall, domain rearrangement involves the displacement of domain I, relative to domains II and III (Fig. 5a/b). One may alternately describe the motion as a relative displacement the II/III superdomain, though structural considerations suggest this process is sterically hindered by the ribosome. 45 In terms of biological function, the balance between EF-Tu dissociation and conformational dynamics is central to aa-tRNA "proofreading," where incorrect tRNA molecules may be rejected, even if they successfully pass the initial association step.</p><p>While the ability of EF-Tu to adopt competing conformational states has been known for decades, 46 the precise physical-chemical relationship between large-scale domain motion and the fidelity of aa-tRNA selection is not fully understood. Recent simulations, 47,48 as well as experiments, 49 have suggested that the timing of EF-Tu release from the ribosome may influence the accuracy of aa-tRNA selection. Experiments have shown that domain rearrangement in EF-Tu and dissociation from the ribosome are kinetically competitive processes. 50 Previous simulations with explicit-solvent models 51 have implicated transient domain separation during this rearrangement. Similarly, simulations with a coarse-grained model predicted that EF-Tu may adopt a disordered intermediate ensemble during domain rearrangement. 45 Together, these studies suggest that the precise level of proofreading exhibited by the cell will be determined by an elaborate interplay of multiple large-scale processes.</p><p>The complicated character of domain rearrangements in EF-Tu makes it an ideal candidate to test the ability of ELViM to identify parallel pathways during biological function.</p><p>Here, we analyze a recently-reported data set in which over 300 spontaneous domain rearrangement events were observed. 45 In contrast to previous studies, where intuitive choices of reaction coordinates had been constructed, 45,51 we find that ELViM can automatically distinguish between the competing pathways of domain rearrangement (Fig. 5c/d). When using ELViM with the internal distances defined by eq. ( 1) and ( 2), there is clear separation of the endpoint configurations (PRE and POST), as well as two distinct transition-state ensembles (TSE1 and TSE2). To measure the differences in the identified ensembles, we next considered the average distance of each residue pair: r ij n , where i and j are the residue indices and n is the ensemble (PRE, TSE1, TSE2, POST). We find that, on average, almost all distances between domain I and the II/III superdomain are larger in TSE1 and TSE2, than in the PRE and POST conformations (not shown). This is consistent with the extended character of the TSEs, where the domains are transiently separated. 45,51 Relative to TSE2, we find that there is a combination of extended and contracted distances in TSE1 (Fig. 5d). Specifically, we find that the distances between the II/III superdomain and some Domain I residues (1, 2, 30-33 and 64-66) are decreased, while nearly all other distances are increased. This implicates relative rotation of the domains, which is consistent with visual inspection of the ensembles (Figure 5e/f). Overall, there is a roughly 180 • rotation of domain I, which is consistent with previous descriptions that implicated extended (TSE1) and disordered (TSE2) ensembles.</p><p>In summary, ELViM is able to automatically detect parallel pathways and competing transition state ensembles, without the need for manual intervention. Rather than having to intuitively design coordinates, which can involve visual inspection of large numbers of transition events, the ELViM method provides a means for efficiently mapping structurallydistinct pathways during elaborate conformational rearrangement in multi-domain, or multicomponent, assemblies.</p><!><p>In this work we presented a method to visualize biological macromolecule energy landscapes.</p><p>The method ELViM is based on the assumption that energy states are highly correlated with its molecular conformations. If one can have a suitable measurement of distances between any two configurations, it is possible to have a proper multidimensional representation of its energy landscape that may then be reduced to a more comprehensible 2D or 3D perspective.</p><p>We have shown that, for the systems analyzed, a suitable metric is based on the inner distances between its atoms.</p><p>While we have focused our analysis on molecular dynamics trajectories, ELViM is applicable to any sampling technique, such as Monte Carlo or Discrete Path Sampling. 52,53 Even though the method has been targeted primarily to proteins, we anticipate this strategy will prove useful for other molecular systems, such as RNA, DNA, or biomolecular assemblies. Further, ELViM is not restricted to any particular level of molecular resolution (C α or all-atom). Rather, the strength of the method relies on the way the distances between any two conformations are measured, where we have found a suitable metric is q-Wolynes (q i,j ), eq.( 1). This metric was initially developed to compare conformational differences between predicted protein structures. In principle, other metrics could be considered for other systems. One strength of ELViM is that it does not require or assume knowledge of an appropriate reaction coordinate. This is particularly beneficial when studying a new system for which minimal quantitative analysis is available.</p><p>Our results suggests a variety of potential applications for this method. For instance, the detailed information of transition state paths, given by the contact analysis and typical structures at TS, can be used to engineer mutations to favor one TS with respect to another. 37 We also envision it will be useful for studying more complex proteins, such as those that contain topological knots, 54,55 or proteins where backtracking mechanisms may occur. 56,57 In closing, we should note that the metric used in this study is based only on structural features, and the mapping of the landscape was not based on kinetic properties. However, there are many avenues by which to further extend this technique, such as including kinetic information to weight the effective distances between conformations, which may enable a more broad thermodynamic and kinetic view of biomolecular dynamics.</p><!><p>The visualization methodology is based in four steps: (I) Generate an ensemble of structures, which in general are obtained through simulations; (II) Calculate the dissimilarity matrix by applying a metric throughout the simulated trajectory; (III) Clustering data, when there is a large number of structures, a data processing procedure is carried out to cluster very similar structures into a single conformation, which represents each cluster; (IV) Perform a multidimensional projection, which transforms the dissimilarity matrix into a 2D projection.</p><p>Each step is discussed in detailed below.</p><!><p>As a first step, one obtains an ensemble of structures. This step is not part of the ELViM itself, in the sense that many techniques may be used to generate the data ensemble. In protein folding investigations, typically one uses a C α , C α + C β , or all-atoms representation. However, ELViM can, in principle, be used in any molecular system; for example it can be used to study energy landscapes of RNA, DNA, macromolecular assemblies and chromatin.</p><p>In this work, we use a structure-based model that represents each protein residue by a single bead placed at the position of its Cα atom. 11 This model defines an experimentally derived structure as the global potential energy minimum. 58 All simulations were performed using the Gromacs (v4.5.4) software package, 59 with force field files generated by the SMOG web server, 60,61 and the native contact map made by CSU algorithm. 62 We performed two sets of simulations, the first one was made to find the folding temperature (T f ). We carried out short simulations (5.0 × 10 8 time steps) in a range of temperature, with T f being calculated using eSBMTools. 63 With T f defined, we carried out long simulations (5.0 × 10 10 time steps) in this temperature. This second set of simulations is performed in order to have the conformational space as thoroughly visited as possible.</p><p>Two proteins were chosen to test the method in the context of folding: src SH3 domain (pdb 1FMK), Protein A (pdb 1BDD); these proteins have been chosen because their folding is well known and they represent basic structural motifs of proteins. 11,36,41 A third protein, the GTP-binding domain of EF-Tu (elongation factor thermo unstable) was also used, since it undergoes a dramatic and large-scale rearrangement during the transition between GTPbound and GDP-bound conformations transition (EFTu-GDP PDB ID: 1DG1, and EFTu-GTP PDB ID: 5YUM).</p><!><p>A critical feature of ELViM is the definition of a metric that coherently and robustly describes the data. For off-lattice systems, such a metric takes into account the protein structure and establishes a distance based in internal conformational differences. This dissimilarity measurement is inspired on replica symmetry breaking parameter of spin glasses, that has been used to constrain clusters of protein states within a given degree of similarity. 32 In the present method, this metric is extended to describe the distance between any two conformations in the entire configuration phase space. The measure of comparison between conformations k and l is given by:</p><p>where r k(l) i,j is the distance between the atom i and j of the conformation k(l). N is the normalization over all non-neighboring pairs of residues and σ i,j = σ 0 |i − j| , being σ 0 = 1 and = 0.15. The dissimilarity between any two structures is described as</p><p>where δ k,l ranges from 0 (identical) to 1 (very different). The dissimilarity measure is calculated between all the structures of a given trajectory, and it yields a n × n triangular matrix,</p><p>where n is the total number of conformations.</p><!><p>This step is necessary if one has a large number of structures. Currently available multidimensional projection procedures are limited to about 50,000 data points, and a well probed phase space may present a much larger number of configurations. In these cases, it is required that similar structures are clustered into single representative configurations prior to performing multidimensional reduction. While this reduces the computational requirements, reducing the number of points does not lead to a lose of information, the number of states associated with each cluster is recorded, and all states of the initial data set are considered in the analysis.</p><p>This procedure is consistent with the other steps of the method, since the same distance criteria between conformations is used throughout the entire method, eq.( 2). Conformations close to each other by less than a cut-off distance δ cut are clustered into a single conformation.</p><p>Starting with a very small δ cut , we group the data so that conformations close to each other (δ k,l < δ cut ) belong to the same cluster. Through iterations δ cut is gradually increased. A condition is set such that if two conformations have a small distance (δ k,l ≤ 0.1) and the difference between distances against all other structures are smaller than ±0.05, this pair of structures is clustered. The clustering is performed until it reaches a stop parameter, which can be a maximum cut-off distance δ M AX cut or a defined number of structures that will be used in the next step. In this way, we preserve all the phase space sampled in the simulation, and at least one representative structure of each region is used to generate the landscape visualization. At the end of this step we have a set of distances between each cluster δ c , which consist of a N × N triangular matrix, where N is the total number of clusters.</p><!><p>The last step in the ELViM method is the reduction of the original multidimensional space, to a low-dimensional (2-3) representation. Multidimensional projection techniques address this question by mapping the data into p-dimensional spaces, with p = {1, 2, 3}, preserving the maximum information about distance or dissimilarity between the data. 64 A multidimensional projection technique can be defined as follows: Let X be a set of N objects in R m , with δ : R m × R m → R a criterion of proximity between two objects in R m , and Y a set of objects R p with p = 1, 2, 3 and d : R p × R p → R a criterion of proximity in R p .</p><p>A multidimensional projection technique can be described as function f , where f : X → Y which aims to make | δ(x i , x j ) − d(x i , x j ) | as close as possible to zero, ∀ x i , x j ∈ X. In our case, this consists of a minimization procedure in which all distances between clusters in the real multidimensional phase space δ k,l correspond to distances d k,l in a 2D phase space.</p><p>These methods are generally known as Multi-Dimensional Scaling methods (MDS). 65,66 It is important to note that ELViM approach is very different from some dimension reduction methods such as PCA, 67,68 in which one seeks component vectors to describe a data. Typically, the variables associated with these vectors are dihedral angles (Φ and Ψ) and specific distances between monomers of atoms.</p><p>MDS procedures usually vary in the way the pairwise distances among a set of N conformations are computed. The classic MDS approach uses a Euclidian distances in the high dimensional space, 69 and Metric MDS methods generalize optimization procedures using different cost functions. 66 Non-linear methods can be broadly classified into two groups: those that provide a mapping (either from the high-dimensional space to the low-dimensional embedding or vice versa), and those that just give a visualization (this is case for the ELViM method). The basic concept of this method is based on proximity data from a dissimilarity matrix (distance matrix). Given a chosen p dimension (usually 2 or 3), the algorithm places each object of the matrix in a position of the p-space dimension, then a minimization is made in order to preserve the distances contained in the original matrix. There are a variety of ways that this procedure is performed, in general, a pre-step is performed in the original matrix, followed by a classic MDS, for example, Isomap 70 uses the Floyd-Warshall algorithm to estimate the geodesic distance matrix, followed by classic MDS to project it into a low dimension; Locally-Linear Embedding (LLE) 71 is very efficient when using sparse matrix because it performs a close neighbor check before performing dimensional reduction.</p><p>The minimization procedure is analogous to the traveling salesman problem, 72,73 in which given a set of distances, one tries to minimize the path to go through all of the cities, which is the cost function. In ELViM, the cost function E is associated with the position of the points in the 2D representation s k,l , and it is given by</p><p>where the sum is carried out over all k and l pairs of conformations.</p><p>Though this step is also an essential part of the ELViM, it can performed using a variety of minimization multidimensional projection techniques. The simplest method of them would be following a Monte Carlo/annealing, but it would be inefficient and time consuming. In the present work, we use the Force-Scheme technique, 64 which is a method based on a standard multidimensional projection technique. It is computationally fast and allows to work with sufficiently large number of conformational states (up to 50,000 points). It offers a balance between precision and computational performance. This iterative minimization technique establishes an approach in which the objects (conformations) are initially positioned, either in a random or heuristic way. Iterations are then performed to introduce perturbations in the positions of each object until the system reaches equilibrium; further details about this technique can be found in our recent work 31 and discussion on system size restrictions can be found in Tejada et. al. 64</p>