Using Dimensionality Reduction to Systematically Expand

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Using dimensionality reduction to systematically expand conformational sampling of intrinsically disordered peptides Oleksandra Kukharenko, Kevin Sawade, Jakob Steuer, and Christine Peter J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00503 • Publication Date (Web): 02 Sep 2016 Downloaded from http://pubs.acs.org on September 13, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Using dimensionality reduction to systematically expand conformational sampling of intrinsically disordered peptides Oleksandra Kukharenko,∗ Kevin Sawade, Jakob Steuer, and Christine Peter∗ Theoretical Chemistry, University of Konstanz, 78547 Konstanz, Germany E-mail: [email protected]; [email protected]

Abstract One of the approaches to improve our ability to characterize biologically important processes and to map out an underlying free energy landscape is to direct MD simulations to explore molecular conformational phase space faster. Intrinsically disordered systems with shallow free energy landscapes of a huge number of metastable minima pose a particular challenge in this regard. Both characterization of the often ill-defined conformational states as well as the assessment of the degree of convergence of phase space exploration are problematic. We have used a multidimensional scaling-like embedding (sketch-map) to describe the energetically accessible regions of phase space for a peptide fragment of the intrinsically disordered protein α-synuclein. Using sketchmap coordinates from a short initial simulation, we guided additional MD simulations to efficiently expand sampling of the conformational space. The sketch-map projections are very well suited to detect rare but possibly functionally relevant events, metastable intermediates and transition states in the vast amount of data.

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Introduction

Intrinsically disordered proteins and peptides (IDPs) are biomolecules that lack well-defined secondary and tertiary structure. The great interest in this class of systems is motivated by an increasing awareness of their biological relevance 1,2 . For example, several proteins that are considered intrinsically unstructured in solution such as α-synuclein, polyQ, τ -protein, and Aβ are involved in neurodegenerative diseases, such as Parkinson’s disease and Alzheimer’s syndrome 3,4 . This is presumably due to their propensity for self-aggregation, fibril formation, but possibly also due to their interaction with biomembranes. 5 Often, these pathophysiologically relevant phenomena are accompanied by conformational transitions where the proteins change from their disordered state into a structurally more defined β-sheet-rich or α-helical state. IDPs pose a challenge for many experimental methods that are traditionally used to study proteins (NMR, 6,7 IR, 8 CD, SAXS, 9 FRET 10 , etc.) and probe well-defined structures or secondary structures. 11,12 For example, NMR spectroscopy is an enormously powerful technique to obtain ensemble-averaged information about the conformations in solution, thus reducing the conformational ambiguity and better define the free-energy landscapes of IDPs, even at amino acid resolution. It is also highly suitable to investigate dynamics. Yet, the characteristic time scales of IDP conformational transitions, the investigation of individual molecules or the molecular mechanisms controlling protein folding upon binding remain a challenge. 7 Here, atomistic molecular dynamics (MD) simulations can be a substantial supplement. 13,14 Assuming ergodicity, MD produces a properly weighted statistical ensemble (e.g. canonical) of conformations and a free energy landscape which can be used to characterize the states and transitions of disordered proteins. 15,16 The ergodicity assumption, however, is problematic in case of IDPs. They are characterized by a comparatively shallow free energy landscape with multiple local minima lacking a dominant global minimum, resulting in a multitude of low-free energy states. In spite of comparatively low-energy barriers, their dynamics for transitions between different metastable states can be slow because of the enormity of con2


Page 2 of 27

Page 3 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


formational space. To obtain a convergent, statistically relevant picture of the dynamics, i.e. an accurate representation of the Boltzmann-weighted distribution of conformations, extensive sampling is required. A number of methods has been developed to enhance the sampling ability of classical MD simulations. In one class of techniques a simulation bias is added in order to force the interesting, rarer events to occur more often and to promote folding/unfolding. Such methods include simulated annealing, 17 replica exchange molecular dynamics (REMD), 18 or replica exchange with solute scaling, 19,20 which will improve the sampling, but not actively drive the system away from already explored regions of phase space. Other methods such as reconnaissance metadynamics, 21 local elevation, 22 NMR-guided metadynamics, 23 etc. 24 explicitly ”push“ the system from already known states. These latter methods typcially require a priory definition of suitable reaction coordinates. The simulation bias forces the value of some particularly interesting collective variables (CVs) to change more rapidly than it otherwise would. While this is more or less straightforward in the case of a well-defined folded structure, the choice of appropriate CVs is not obvious in the case of partially ordered systems. A different category of methods such as Markov State Models (MSM), 25 directional milestoning, 26 forward-flux sampling, 27 known as state-based methods, combine information from shorter simulations between various locations in phase space to obtain global information such as reaction paths and rates, instead of relying on the statistical analysis of long trajectories. As such, they may be more advantageous in systems with many metastable minima and multiple paths. 24 However, they require the initial set of conformations to cover a wide area of the energy landscape. Especially for shallow free-energy landscapes with a vast amount of metastable states this can be problematic. For such systems, a need of sampling algorithms arises, that neither rely on prior knowledge of the states the system visits nor on predefined reaction coordinates. In this paper, we propose a sampling scheme which uses a multidimensional scaling-like (MDS) embedding (sketch-map), 28,29 to characterize and effectively explore the conforma-

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 4 of 27

Page 5 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2 2.1

Methods System and Simulation Details

To illustrate the method that will be explained in Section 2.3 we chose a 22-residue long peptide fragment (42-63) of α-synuclein, which, as suggested in ref. 33, includes a fouramino acid linker important for fibril formation. Its sequence in one-letter code is SKTKEGVVHGVATVAEKTKEQV. Since α-synuclein is highly interface active, with conformational transitions that are driven by interface interactions, the peptide was placed at a water/vacuum interface. The presence of the interface further slows conformational dynamics, which pose additional sampling challenges and makes the chosen system and setup an ideal test case for our algorithm. All simulations were carried out on an atomistic level using the GROMACS 4.6.5 program package 34 , SPC water 35,36 and the GROMOS 54A7 force-field 37 . Electrostatic interactions were handled by the particle mesh Ewald method (PME) 38,39 with a real-space cut-off distance of 1.0 nm and a Van der Waals cut-off distance of 1.4 nm. A standard molecular dynamics simulation was performed using the leapfrog integrator with 2 fs integration steps and a velocity-rescaling algorithm 40 for the temperature coupling at 300 K. No pressure coupling algorithm was introduced to ensure a stable interface. The system was prepared with the peptide immersed in the water phase close to the surface. The peptide moved to the interface and stayed there during the simulation – with only very brief intervals of (partial) immersion. Simulations were performed on 16 processing cores/32 threads, with 12800 CPUh per µs simulation time.

2.2

Sketch-map

Dimensionality reduction is a key element in identifying the regions of conformational phase space the system has visited and in subsequently grouping similar structural motives in the trajectories (clustering) and vice versa in determining sparse (low-probability) regions. Direct clustering of high-dimensional data usually faces a number of challenges: a concept of 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 27

distance becomes less precise as the number of dimensions grows, meaning the notion of nearest and farthest point becomes meaningless and the data become sparse (known as ”the curse of dimensionality“); not all dimensions are uncorrelated (i.e. they do not add any new information to the description); and it’s harder to identify similarities in systems with many characteristics (dimensions). Therefore, for clustering of high-dimensional data usually some dimensionality reduction is used, either by projection the data to low-dimensional space or by clustering of subspaces of high-dimensional space (simple examples are clustering based on RMSD, binning the Ramachandran plot for small peptides, etc). Here to reduce the number of dimensions we used the sketch-map algorithm, 29 which can project highly nonlinear conformation space. The high-dimensional coordinates that accurately characterize the conformations can be internal distances, Ramachandran angles, or other sets of collective variables. The algorithm is based on metric multidimensional scaling and its advantage is that it focuses on an intermediate range of similarities/dissimilarities. It disregards information that corresponds to thermal fluctuations around (meta)stable structures as well as exact reproduction of large values of dissimilarities caused by the high-dimensionality of space. It generates a set of low-dimensional projections from a set of high-dimensional landmark points (a small subset of the simulated data) by minimizing the function of dissimilarities of distances between these points

χ2 =

X

[s(Rij , A, B) − s(rij , a, b)]2 ,

(1)

i6=j

where Rij and rij are dissimilarities (e.g. Euclidean distances) between high-dimensional landmark points Xi and Xj and their low-dimensional projections xi and xj , respectively; s(r, a, b) is the following sigmoid function s(r, a, b) = 1 − (1 + (2a/b − 1)(r/σ)a )−b/a ,

6


(2)

Page 7 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


which transforms all distances to values between zero and one. Selection of sigmoid function parameters σ, A, B, a, b (which should save only important features of a phase space) is done by analyzing the distance distribution between high-dimensional points. Using the projection of the initial landmark projections of all other points from highdimensional space, X, are found by iteratively minimizing a strongly non-linear fitness function

σ 2 (x) =

N X

[s(Ri (X), A, B) − s(ri (x), a, b)]2 → min,

(3)

i=1

where Ri (X) is the dissimilarity between X and i th landmark point and ri (x) is the distance between point x and the projection of the i th landmark.

2.3

Expansion

The sketch-map projection generates dense basins from structurally close conformations. Therefore, rarely visited distinctly different conformations as well as conformations that lie on transition paths between basins can be easily identified. By selecting starting conformations from the sparsely populated regions we can increase the probability that the system explores so-far unvisited regions of phase space. It was suggested, 29 that the sketch-map algorithm will generate meaningful coordinates even if the already mapped trajectory does not visit all the energetically accessible parts of conformational space, so all the additional simulations can be projected using the initially obtained and projected small subset of representative points (landmarks). The landmark points do not have to be necessarily present in all simulations. They can be generated separately, selected from initial simulations and the set can also be extended during the broadening of conformational space. (Note: If selected landmarks poorly represent new conformational space then either additional landmarks from new projections or totally new landmarks can be selected from all available trajectories). The expansion scheme can be summarized in the following algorithm:

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(1) Run initial MD simulation(s) from the starting conformation(s). (2) Select landmark points from the simulation(s) and perform dimensionality reduction using the sketch-map algorithm. Project all other points from the simulation and obtain a low-dimensional representation of the free energy surface as a function of sketch-map coordinates. (3) Select N random points from the low-density regions of the obtained projection. (4) Run new MD simulations from the selected conformations with new initial velocities and a short (a few picoseconds) position restraining of the peptide. (5) Make projections of new simulations using landmark points from step 2. Repeat step 2 with additional landmarks if sampling appears insufficient. (6) Repeat steps 3-5 until the sampled phase space converges: new states are not obtained. Steps 4 and 5 can be run simultaneously (projections can be calculated on the fly). There is a possibility to use certain stopping criteria at step 4: a simulation could be terminated in two cases: (a) if it falls back to already explored states or (b) it has formed a new state and stays there for a predefined amount of time. Then the method is divided into two stages: exploring and refining (with help of other enhance sampling techniques or without, as shown in ref. 30–32). For the present case we have decided to run many parallel, comparatively long simulations without introducing any bias. This allows to directly use the obtained data for further analysis, e.g. to extract kinetic properties afterwards using methods such as MSM and to estimate free energies for the states of interests.

8


Page 8 of 27

Page 9 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3

Results

3.1

Applying the algorithm

As an initial step we performed a 1 µs simulation of the 22-residue α-synuclein fragment starting from a helical conformation. For analysis we selected frames every 10 ps (in total 1 × 105 snapshots). Sketch-map coordinates were generated by selecting 2000 landmark points from the trajectory using a farthest point strategy. The full set of pseudo-torsion angels between Cα atoms were used as the high-dimensional collective variables – resulting in 20 dimensions. Euclidean distances (taking into account the periodicity of torsional angle values) were used as dissimilarity measure in high and low dimensional space. The optimal two-dimensional projection was obtained with the following parameters for the sigmoid function σ = 5, A = 20, B = 6, a = 2, b = 6 (see Eq. 2). The obtained projection is shown in Figure 1a. The figure nicely demonstrates one difficulty when one simulates IPDs. How does one assess convergence, i.e. exhaustive/sufficient exploration of phase space for a system which does not possess well-defined structures but which is not entirely structureless either? How does one know that one has found all relevant states? For the (patho)physiological function of these systems it may be essential to capture certain low-populated states of partially formed structures. Some structures may be relevant in conformational-selection scenarios when the systems undergo transitions to β-sheets in amyloid aggregates or to α-helices when interacting with membranes. The method suggested here can help to address these problems. Moreover, analysis of some clusters grouped together by sketch-map due to structural similarity (see snapshots in Figure 1a) suggests that neither a conventional clustering analysis based on atom-positional root-mean-square deviations nor a secondary-structure based analysis would easily have identified all those states. From the initial simulation shown in Fig. 1a we started the expansion. From the lowdensity regions we randomly selected new 100 starting points (Fig. 1b) for the expansion in step 4. The new simulations amounted to a total of 10 µs, and again timeframes every

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

10 ps were selected for the projection (i.e. a total of 11 × 105 analyzed snapshots). The same landmarks and parameters of the sigmoid function were used. After the analysis of the obtained projection step 2 was repeated: new 2000 landmarks were selected from all obtained trajectories and a new dimensionality reduction was conducted (with the same sigmoid function parameters). A further extension was prepared with 100 new simulations and an aggregated simulation time of 3 µs. The collection of landmark points was extended to 2500. The final projection as well as the evolution of the sampled configurational phase space can be seen in Figure 2. The projection in the inset nicely shows that the system has effectively explored new regions of conformational space and sampled new structures that had not been visited during the initial simulation. Note, that due to the bias in selecting conformations from sparse regions for the extension with the main aim to explore the accessible configurational phase space, this final projection is in fact not a properly weighted depiction of the free energy landscape of the system. Nevertheless it is useful, to display the sampled landscape in terms of Boltzmann inverted probabilities to identify the centers of local free energy minima and to define structures or states for subsequent analysis. The figure also shows representative structures of local free energy minima (with numbering indicating their population). Minima in the free energy surface were automatically detected by binning the two-dimensional projection to rectangles (cells) and selecting iteratively local maxima (most populated cells, separated by change in the density gradient). Vice versa, bins that have only one point in it and are surrounded by empty bins are called low-density or sparse regions. From the latter new starting configurations for step 3 are selected. During the expansion, the size of the projection usually increases, but the size of the bins is kept constant, only their number increases. It is possible to use different definitions for dense/sparse regions, e.g. introducing some cut-off distance for points and counting neighbors within this cut-off. For minima showed in Figure 2 we have analyzed which simulation trajectories contributed to them and whether minima have been repeatedly visited by several separate

10


Page 10 of 27

Page 11 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


trajectories. This allows a first, cautious assessment of the convergence of the approach. To illustrate several scenarios that can occur during the expansion, the trajectories which contributed to minima 1, 3, and 5 and the relative time they spent there is shown in Figure 3. The deepest (first) minimum was present at the initial simulation and was extensively populated during the expansion stage by five more trajectories (which got ”trapped“ in this minimum). It would have been possible to use the stopping criteria and prevent the simulations to fall back into this minimum, thus making the exploration stage more efficient, but since we plan to subsequently use the data of this system for kinetic analysis, we decided not to use any stopping criteria here. The third minimum was populated by the initial simulation as well as by several new ones, but in this case simulations were able to ”escape“ this metastable state. There are also minima, which are populated only by new simulations (e.g. the fifth deepest minimum, see Fig. 3).

3.2

Zoom in and automated event detection

The free energy landscape of IDPs appears qualitatively inverted with respect to that of structured proteins. 23 Meaning that the deepest minima will be populated by unstructured motifs, while partially folded conformations may appear as transitional states with relatively similar free energies and rapidly interconnect among each other. Detection of these events is important, since they may be responsible for signaling and regulatory functions of proteins and possibly important in view of mechanistic studies of environment-induced conformations. 41 The advantage of using sketch-map coordinates for objects such as IDPs is the possibility to use other descriptors (e.g. time, temperature of a replica, some energetic quantity, etc.) to pick small (but highly relevant details), which are easily overlooked. By automatically analyzing distance distributions in the short trajectories (or trajectory fragments) it is possible to detect metastable states and transition states. Fig. 4a nicely illustrates this by showing a short trajectory fragment in which a conforma11




Page 12 of 27

Page 13 of 27




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tional transition occurs, as can be seen both in the (zoomed in) sketch-map projection and the conventional secondary structure plot. In such a case the distribution of all (low-dim.) distances between pairs of conformations is bimodal, and the integral of the distribution exhibits distinct steps. Such distinct steps in the integral over the distance distribution indicate that a transition between states has happened. Transition points can be identified as those points where the difference in the distances to the previous and the subsequent points (as function of time) is significant. Fig. 4b presents a second example of a successful detection of an important event that might easily have been missed among a large amount of data. We have identified the formation of a full α-helix, which has happened only once in 14 µs of simulation time. In addition, by using not only density information, but also the colouring based on the secondary structure as a additional description for each point (colouring scheme is the same as in Fig. 1a), one can identify regions of interest, even if they are poorly populated, and automatically select new simulations which contribute to those regions (although they were started in different locations of phase space). It is straightforward to use the basins of the sketch-map and run transition path sampling simulations if one is interested in this specific event for the current IDP.

3.3

Pushing the algorithm

So far, the algorithm has been shown to work and to harness sketch-map to automatically assess and improve the sampling of conformational phase space of IDPs and to detect events even if these are only spurious folding events into partially-formed, short-lived structures. However, the simulation data (a total of 14 µs of simulation time) presented so far do not yet conclusively show that the expansion algorithm – if brutally used for expansion only – is indeed significantly more efficient than brute-force MD. In the following example, the same α-synuclein peptide fragment had been prepared in a β-sheet starting conformation, immersed in water and moved to a water/vacuum interface. 14


Page 14 of 27

Page 15 of 27




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

After 100 ns of simulation time it turns out, that this starting conformation leads to a kinetically trapped strand-loop-strand structure at the interface with frequently forming and loosening β-sheets. In 1 µs MD simulation time, no significant conformational change could be observed (see Fig. 5a). We used the expansion algorithm with stopping criteria to compare, how much of the conformational space could be recovered by the β-sheet simulation at the same simulation time, and whether the same metastable states as in the previous 14 µs simulation (started from an α-helix) are visited during the expansion. The algorithm was applied in the following way: The first 100 ns of the simulation (bulk water phase) were fully excluded from the analysis. The following 100 ns were considered as an initial simulation (i.e. as step 1 of the algorithm). This means all new simulations were conducted at the water/vacuum interface. After projecting the initial simulation, 500 landmarks were selected. During the expansion steps this number was increase to 1500 points. Sigmoid function parameters were changed as the distance distribution profile was changing while new conformations were explored. At step 3, 20 random frames were chosen from sparsely populated regions (they are defined by binning the projection as explained in section 3.1) and new simulations were started. When the new simulations reached 5 ns of simulation time, stopping criteria were applied to each of them with the following termination conditions: if a simulation falls to an already explored state (here we call a state any number of bins having more than a fixed amount of points) and stays there for 1 ns or if the simulation forms a new state, but does not exit it during 0.2 ns, it should be terminated; and a new simulation had to be started from a new random point selected from sparse regions of the joint projection of the initial simulation and all expansion runs. If the stopping conditions were not fulfilled, the simulation was continued and checked again each 5 ns. We stopped the expansion, when we reached a total 900 ns of simulation time. The comparison of the projections of an initial 1 µs simulation (bulk phase excluded) with 1 µs aggregated simulations obtained using the algorithm (composed of 100 ns initial simulation and the 900 ns expansion) are shown in Figure 5. It is evident, that during the

16


Page 16 of 27

Page 17 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


same simulation time we were able to cover a much larger area of conformational space. By comparing this simulation with the one described in section 3.1, we see that in this comparatively short simulation time, we were able to ”visit“ most of the deepest minima, detected in 14 µs of simulation (Fig. 5c). We would like to point out that Figure 5c could be misleading with respect to the continuity of trajectories. The blue points are a projection of a new simulation into a sketch-map of a previous one, i.e. the blue trajectories have not been well-represented by the landmarks used to span the space. Therefore this Figure cannot necessarily be used as an indication that sketch-map disrupts trajectories. It is only meant as proof, that the drastic expansion method had succeeded at ”blasting“ the system out of the deadlock structure (hairpin) at the interface, therefore we did not re-sketch-map the two ensembles together (which would have avoided the strong impression of discontinuity, and allow to observe traces of small conformational transitions in the sketch-map projection). In general projections to low-dimensional space are not necessary continuous even though by construction the trajectories generated by the dynamics are continuous and the expansion steps start in regions of conformation space that had been visited by a previous trajectory. It should be noted, that due to its focus on preserving intermediate distances the sketch-map algorithm has been found to be well-suited to represent conformational transitions and continuity in phase space 29 .

4

Discussion

The previous sections show that sketch-map is excellently suited to reduce the dimensionality of conformational space and to project and group data in the case of intrinsically disordered systems with shallow free-energy landscapes. Although other valuable methods for dimensionality reduction exist, we expect sketch-map to be particularly suitable for our expansion scheme compared to other non-iterative linear methods, such as principal component analysis 42 or time-lagged independent component analysis (TICA) 43 or other non-linear methods

17




Page 18 of 27

Page 19 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(locally linear embedding, 44 Isomap, 45 and diffusion maps 46,47 ). For example, TICA uses the kinetic information of the system for extracting the slowest collective motion. This is a good choice for systems with a small number of well defined metastable states, but all fast transitions will be inseparable, at least at the small number of dimensions. For the present (IDP) system we found that to obtain a reasonable projection using TICA one needs more than 10 dimensions. In low-dimensional space with significantly more dimensions than two, identifying truly sparse regions with a limited amount of data points is difficult. Diffusion map was successfully used for the exploration of conformational space for several systems with a limited number of states/metastable intermediates, 30,31,48 but it requires to calculate the distance matrix (and its eigenvalues) between all the frames in all trajectories at each extension step to perform dimensionality reduction. This might be challenging in terms of computational time with the growth of trajectory length. In ref. 30,31 the locally scaled diffusion map was used to select starting conformations for new simulations in the exploring stage based on diffusion coordinates calculated as eigenfunctions of the kernel, which is derived from the root-mean-square displacement (RMSD) between all conformations and ”local scale“. Determining the last parameter requires a Boltzmann-distributed set of molecular configurations and the performance of MDS for each conformation within some neighborhoods. 49 In ref. 32 a MSM was used to guide the exploration stage: At each step of clustering RMSD calculations were conducted to define Markov states and a new starting point using the mean first passage time of obtained MSMs. For IDP system with an extremely high number of metastable states, using RMSD as the main dissimilarity criterion is problematic, the here presented MDS-based algorithm may have advantages. Applying abovementioned methods may also be challenging, because there may not be a significant difference in the eigenvalues, meaning there is no timescale separation which is important for MSM or big structural difference determined by diffusion map. Additionally sketch-map uses only a small number of reference structures (landmarks) to produce good out-of-sample embedding. It’s worth to mention, that these coordinates

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

can be used as collective variables to construct bias potentials to enhance sampling 50 . The proposed scheme can be embedded in automatic adaptive sampling packages, e.g. HTMD 51 . At present we have used MDS-based characterization only for expansion and exploration. We envision this method as a very first step in a more quantitative process towards a proper free energy landscape. Here it merely explores the conformational ”landscape“ (heuristically and qualitatively). During this process, the landscape is dynamically clustered into states, and transitions are recorded. Thus, the expansion scheme is rather a search scheme than a sampling scheme, but it can be extended to become a sampling scheme by combining with TPS-like or forward-flux methods 27 with sketch-map coordinates as the reaction coordinates to obtain more quantitative estimates of energy differences between states determined by our expansion scheme. A more rigorous use of the method to generate free energy landscapes, using reweighting algorithms such as TRAM, 52 is planned. 48,52–54 With the method proposed here, we were trying to address a particularly formidable challenge in the simulation of IDPs caused by their high conformational heterogeneity. We need to stress that given the limitation of atomistic forcefields, comparison with experimental data is crucial to asses the physical relevance of the simulated data. Here we have refrained from interpreting the results with regard to the particular system and the forcefield used. Previous force field comparison studies have demonstrated a high sensitivity of IDP conformational ensembles to differences between the force fields 55,56 . There is a substantial interest in obtaining accurate ensembles of IDPs using all-atom simulations. We believe, that our algorithm can be also very efficiently applied to a qualitative comparison of conformational phase space sampling of IDPs in different forcefields. 55,57 Furthermore, the data of the simulated ensembles can be then used for a quantitative comparison with experimental methods that give access to for example distance distributions such as NMR or EPR spectroscopy.

20


Page 20 of 27

Page 21 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5

Conclusions

We proposed an expansion scheme to efficiently explore conformational space, which uses the sketch-map dimensionality reduction algorithm. An advantage compared to many other methods is that one neither needs to have a general knowledge of the states nor a reaction coordinate. The algorithm is best suited, if one knows rather little of the conformations the system might eventually adopt, and if it is very complex, i.e. while the sampling is expanding the picture of the conformational landscape of the system changes drastically. It allows to project newly sampled areas to the known landscape to assess the expansion. At the same time, since the algorithm allows resketching whenever necessary, the two coordinates that the ensemble is projected into change adaptively, resulting in an optimal redistribution to evaluate sparsity and find new points from where to expand. This distinguishes the method from other projection methods, where the a priory decision on the coordinates of the projection is a limitation (especially if one does not know where the system will go eventually). Here, the coordinates are optimally chosen depending on the current state of sampling. The algorithm is particularly suited for disordered systems with many metastable, partially structured states. It was successfully applied to a 22-residues peptide fragment of α-synuclein at a water/vacuum interface. This is a challenging system because of the vast number of conformational states and extremely long lived metastable states populated by unstructured or partly structured motifs stabilized by the interface. At the same time, folded conformations, which are relevant for comparison with experimental measurements are rare and unstable and cannot be easily distinguished by most methods relying on the dynamics of the system. It can effectively distinguish between different conformational motifs and separate them in two-dimensional space using values from high-dimensional conformational descriptors, such as the torsion angels between Cα atoms. Projections of new data from the ongoing simulation can be done on the fly. The exploration of the accessible phase space can be very efficiently guided. Obtained data can be directly used for further analysis such 21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

as MSM, and does not require further refinement or an unbiasing stage.

6

Acknowledgements

We would like to thank I. Kevrekidis for very helpful discussions. This work was supported by a fellowship of the Zukunftskolleg of the University of Konstanz to OK and by the SFB969 funded from the German Science Foundation. Simulations were performed on the computational resource bwUniCluster funded by the Ministry of Science, Research and the Arts Baden-W¨ urttemberg and the Universities of the State of Baden-W¨ urttemberg, Germany, within the framework program bwHPC.

References (1) Tompa, P. Trends Biochem. Sci. 2002, 27, 527 – 533. (2) Dyson, H. J.; Wright, P. E. Nat. Rev. Mol. Cell Biol. 2005, 6, 197–208. (3) Uversky, V. N.; Oldfield, C. J.; Dunker, A. K. Annu. Rev. Biophys. 2008, 37, 215–246. (4) Goedert, M. Science 2015, 349, 601. (5) Uversky, V. N.; Dunker, A. K. In Protein and Peptide Folding, Misfolding, and NonFolding; Schweitzer-Stenner, R., Ed.; John Wiley and Sons, Inc.: Hoboken, N. J., 2012; pp 1–54. (6) Jensen, M. R.; Ruigrok, R. W.; Blackledge, M. Curr. Opin. Struct. Biol. 2013, 23, 426 – 435. (7) Jensen, M. R.; Zweckstetter, M.; rong Huang, J.; Blackledge, M. Chem. Rev. 2014, 114, 6632–6660.

22


Page 22 of 27

Page 23 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(8) Natalello, A.; Ami, D.; Doglia, S. M. In Intrinsically Disordered Protein Analysis: Volume 1, Methods and Experimental Tools; Uversky, N. V., Dunker, K. A., Eds.; Humana Press: Totowa, NJ, 2012; Chapter Fourier Transform Infrared Spectroscopy of Intrinsically Disordered Proteins: Measurement Procedures and Data Analyses, pp 229–244. (9) Bernado, P.; Svergun, D. I. Mol. BioSyst. 2012, 8, 151–167. (10) Hofmann, H.; Soranno, A.; Borgia, A.; Gast, K.; Nettels, D.; Schuler, B. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 16155–16160. (11) Uversky, V., Longhi, S., Eds. Instrumental Analysis of Intrinsically Disordered Proteins: Assessing Structure and Conformation; John Wiley and Sons, Inc.: Hoboken, N. J., 2010. (12) Uversky, N. V., Dunker, K. A., Eds. Intrinsically Disordered Protein Analysis: Volume 1, Methods and Experimental Tools; Humana Press: Totowa, NJ, 2012. (13) Knott, M.; Best, R. B. PLoS Comput. Biol. 2012, 8, 1–10. (14) Burger, V. M.; Gurry, T.; Stultz, C. M. Polymers 2014, 6, 2684. (15) Stanley, N.; Esteban-Mart´ın, S.; Fabritiis, G. D. Prog. Biophys. Mol. Biol. 2015, 119, 47 – 52. (16) Qiao, Q.; Bowman, G. R.; Huang, X. J. Am. Chem. Soc. 2013, 135, 16092–16101. (17) Bernardi, R. C.; Melo, M. C.; Schulten, K. Biochim. Biophys. Acta, Gen. Subj. 2015, 1850, 872 – 877. (18) Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 1999, 314, 141 – 151. (19) Liu, P.; Kim, B.; Friesner, R. A.; Berne, B. J. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 13749–13754. 23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(20) Wang, L.; Friesner, R. A.; Berne, B. J. J. Phys. Chem. B 2011, 115, 9431–9438. (21) Tribello, G. A.; Ceriotti, M.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 17509–17514. (22) Huber, T.; Torda, A. E.; Gunsteren, W. F. J. Comput. Aided Mol. Des. 1994, 8, 695–708. (23) Granata, D.; Baftizadeh, F.; Habchi, J.; Galvagnion, C.; De Simone, A.; Camilloni, C.; Laio, A.; Vendruscolo, M. Sci. Rep. 2015, 5, 15449. (24) Rohrdanz, M. A.; Zheng, W.; Clementi, C. Annu. Rev. Phys. Chem. 2013, 64, 295–316. (25) Bowman, G. R., Pande, V. S., Noé, F., Eds. An Introduction to Markov State Models and Their Application to Long Timescale Molecular Simulation; Springer Netherlands: Dordrecht, 2014; Vol. 797. (26) Májek, P.; Ron, E. J. Chem. Theory Comput. 2010, 6, 1805–1817. (27) Allen, R. J.; Valeriani, C.; ten Wolde, P. R. J. Phys.: Condens. Matter 2009, 21, 463102. (28) Ceriotti, M.; Tribello, G. A.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 13023–13028. (29) Ceriotti, M.; Tribello, G. A.; Parrinello, M. J. Chem. Theory Comput. 2013, 9, 1521– 1532. (30) Zheng, W.; Rohrdanz, M. A.; Clementi, C. J. Phys. Chem. B 2013, 117, 12769–12776. (31) Preto, J.; Clementi, C. Phys. Chem. Chem. Phys. 2014, 16, 19181–19191. (32) Zhou, T.; Caflisch, A. J. Chem. Theory Comput. 2012, 8, 2134–2140. (33) Shvadchak, V. V.; Subramaniam, V. Biochemistry 2014, 53, 279–281. 24


Page 24 of 27

Page 25 of 27

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(34) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435–447. (35) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Intermol. Forces 1981, 331–342. (36) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269– 6271. (37) Schmid, N.; Eichenberger, A. P.; Choutko, A.; Riniker, S.; Winger, M.; Mark, A. E.; Gunsteren, W. F. Eur. Biophys. J. 2011, 40, 843–856. (38) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089–10092. (39) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577–8593. (40) Bussi, G.; Donadio, D.; Parrinello, M. J. Chem. Phys. 2007, 126, 014101. (41) Dalgicdir, C.; Globisch, C.; Peter, C.; Sayar, M. PLoS Comput. Biol. 2015, 11, 1–24. (42) Pearson, K. Philos. Mag. 1901, 2, 559–572. (43) Pérez-Hernández, G.; Paul, F.; Giorgino, T.; De Fabritiis, G.; Noé, F. J. Chem. Phys. 2013, 139, 015102. (44) Roweis, S. T.; Saul, L. K. Science 2000, 290, 2323–2326. (45) Tenenbaum, J. B.; Silva, V. d.; Langford, J. C. Science 2000, 290, 2319–2323. (46) Nadler, B.; Lafon, S.; Coifman, R. R.; Kevrekidis, I. G. Appl. Comput. Harmon. Anal. 2006, 21, 113–127. (47) Kim, S. B.; Dsilva, C. J.; Kevrekidis, I. G.; Debenedetti, P. G. J. Chem. Phys. 2015, 142, 085101. 25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(48) Das, P.; Frewen, T. A.; Kevrekidis, I. G.; Clementi, C. In Coping with Complexity: Model Reduction and Data Analysis; Gorban, N. A., Roose, D., Eds.; Springer Berlin Heidelberg: Berlin, Heidelberg, 2011; Chapter Think Globally, Move Locally: Coarse Graining of Effective Free Energy Surfaces, pp 113–131. (49) Rohrdanz, M. A.; Zheng, W.; Maggioni, M.; Clementi, C. J. Chem. Phys. 2011, 134, 124116. (50) Tribello, G. A.; Ceriotti, M.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 5196–5201. (51) Doerr, S.; Harvey, M. J.; Noé, F.; De Fabritiis, G. J. Chem. Theory Comput. 2016, 12, 1845–1852. (52) Wu, H.; Paul, F.; Wehmeyer, C.; Noé, F. Proc. Natl. Acad. Sci. U. S. A. 2016, 113, E3221–E3230. (53) Hummer, G.; Kevrekidis, I. G. J. Chem. Phys. 2003, 118, 10762–10773. (54) Mey, A. S. J. S.; Wu, H.; Noé, F. Phys. Rev. X 2014, 4, 041018. (55) Henriques, J.; Cragnell, C.; Skepö, M. J. Chem. Theory Comput. 2015, 11, 3420–3431. (56) Rauscher, S.; Gapsys, V.; Gajda, M. J.; Zweckstetter, M.; de Groot, B. L.; Grubmller, H. J. Chem. Theory Comput. 2015, 11, 5513–5524. (57) Best, R. B.; Mittal, J. Proteins: Struct., Funct., Bioinf. 2011, 79, 1318–1328.

26


Page 26 of 27

Page 27 of 27



Using Dimensionality Reduction to Systematically Expand

Recommend Documents