FULL PAPER
WWW.C-CHEM.ORG
Protein Folding Pathways Extracted by OFLOOD: Outlier FLOODing Method Ryuhei Harada,*[a,b] Tomotake Nakamura,[c] Yu Takano,[b,d] and Yasuteru Shigeta[a,b] The Outlier FLOODing method (OFLOOD) is proposed as an efficient conformational sampling method to extract biologically rare events such as protein folding. In OFLOOD, sparse distributions (outliers in the conformational space) were regarded as relevant states for the transitions. Then, the transitions were enhanced through conformational resampling from the outliers. This evidence indicates that the conformational resampling of the sparse distributions might increase chances for promoting the transitions from the outliers to other meta-stable states, which resembles a conformational flooding from the outliers to the neighboring clusters.
OFLOOD consists of (i) detections of outliers from conformational distributions and (ii) conformational resampling from the outliers by molecular dynamics (MD) simulations. Cycles of (i) and (ii) are simply repeated. As demonstrations, OFLOOD was applied to folding of Chignolin and HP35. In both cases, OFLOOD automatically extracted folding pathways from unfolded structures with ns-order computational costs, although ms-order canonical MD failed to extract them. C 2014 Wiley Periodicals, Inc. V
Introduction
used for many targets, such as spin glasses[5] and folding of small peptides.[6,7] Those with the non-Boltzmann weight, such as Multicanonical Monte Carlo[8] and Multicanonical MD simulations,[9,10] were also developed and are commonly used in biophysics for the enhanced conformational sampling of biomolecules. As another powerful structural sampling method, meta-dynamics[11] provides FELs projected onto a subspace spanned by a set of reaction coordinates (RCs) using historydependent repulsive biasing potentials. As strategies for searching the transition pathways, transition path sampling[12] and the string method[13] have rigorous formulations to obtain conformational transition path ensembles. The above conformational sampling methods have been well established and widely used in several biological systems. However, there might exist an issue to be reconsidered when applying these methods to large systems. Generally, with the increase in the vast number of degrees of freedom (DOF), the conformational spaces to be swept become complicated due to high dimensionalities.
Rare events in biological systems are strongly related to their functions, such as ligand binding or allosteric regulation, and so on. In most biological processes, the rare events are observed as structural transitions on biological functions, and thus their extractions are indispensable for understanding biological relevancy. Molecular dynamics (MD) simulation is one of the primary tools to assess biologically important reactions, as a time-series of trajectories at atomic resolution. However, the time scale of biological functions might exceed the accessible time scale of conventional molecular dynamics (CMD) simulations. Recently, a massive supercomputer designed by D. E. Shaw Group has simulated the microsecond-order folding simulations of 12 structurally diverse proteins and unveiled a series of folding mechanisms at atomic levels.[1] From the point of view of computational costs, these long-time CMD simulations seem expensive. In addition, the biologically rare events such as protein folding tend to be observed as stochastic processes. Therefore, it does not ensure that the long-time CMD simulations surely extract the rare events due to the uncertainties coming from the stochastic processes. From the point of the view of free energy landscape (FEL) calculations, CMD simulations with the Boltzmann weight might fail to sample the broad conformational spaces because of trapping into local minima separated by high free energy barriers. To tackle the insufficiency of the structural sampling, developments of methods for promoting the structural transitions are also indispensable for extracting the biologically rare events. Herein, there are two fundamental issues to solve in this field. One is how to find the global and local minima, which are biologically relevant. The other is how to draw FELs. As methods for investigating broad FELs, generalized ensemble algorithms with the Boltzmann weight, such as the parallel tempering[2] and replica exchange method (REM),[3,4] have been widely
DOI: 10.1002/jcc.23773
[a] R. Harada, Y. Shigeta Division of Life Science, Center for Computational Sciences, University of Tsukuba, Tennodai, Tsukuba, Ibaraki 305-8577, Japan E-mail: [email protected] or [email protected] [b] R. Harada, Y. Takano, Y. Shigeta JST-CREST, Kawaguchi, Saitama 332-0012, Japan [c] T. Nakamura Fujitsu Limited, Numazu, Shizuoka 411-0301, Japan [d] Y. Takano Research Center for State-of-the-Art Functional Protein Analysis, Institute for Protein Research, Osaka University, Suita, Osaka 565-0871, Japan Contract grant sponsor: JSPS KAKENHI (Young Scientists (B)); Contract grant number: 26840056; Contract grant sponsor: Scientific Research on Innovative Areas; Contract grant sponsor: MEXT; Contract grant numbers: 26107004 and 26105012; Contract grant sponsor: Computational Materials Science Initiative (CMSI; Consortium) C 2014 Wiley Periodicals, Inc. V
Journal of Computational Chemistry 2015, 36, 97–102
97
FULL PAPER
WWW.C-CHEM.ORG
In treating complicated biological systems, the key point is how to treat the high dimensionalities of the conformational spaces due to increasing of DOF. To describe states of biological systems, a set of RCs should be selected and highdimensional spaces are defined by RCs. Here, to analyze the high-dimensional distributions, clustering techniques might be powerful, as they directly and systematically treat the high dimensionalities. In the clustering processes, stable states of a given protein are regarded as dense distributions in the conformational spaces, that is, clusters. In contrast, sparse distributions are detected as outliers that do not belong to the clusters. These outliers might have potentials to transit to neighboring clusters with high probabilities. Therefore, in this study, we regard the outliers as relevant states that existing in transitional regions. Furthermore, we show that the conformational resampling from the outliers via short-time MD simulations tends to promote the structural transitions of proteins, which efficiently leads to extractions of biologically rare events.
Methods Based on the detections of outliers and their conformational resampling, we propose an efficient conformational sampling method to extract the biologically rare events of proteins, called Outlier FLOODing method (OFLOOD). In OFLOOD, a set of cycle consisting of the following schemes is simply repeated: (i) detections of the outliers from conformational distributions by clustering of MD trajectories and (ii) conformational resampling from the detected outliers through shorttime MD simulations. Here, the detections of the outliers and clustering of MD trajectories are performed using a newly developed clustering algorithm, FlexDice.[14–17] The conformational resampling is performed by restarting short-time MD simulations from the outliers via regenerations of initial velocities. In several past studies,[18–23] restarting MD simulations via regenerations of initial velocities might enhance the conformational sampling of proteins. In every cycle, outliers are updated based on cumulative MD trajectories obtained from the past conformational resampling. These cycles are continued until the conformational search sufficiently converges, which means that the convergences of cumulative distributions projected onto each RC are checked for the termination of OFLOOD. Figure 1 depicts the physical meaning of tracing the outliers in OFLOOD. As shown in Figure 1, the conformational resampling from the outliers might expand distributions through transitions to the neighboring clusters, which leads to extractions of biologically rare events of proteins. The strategy of OFLOOD based on the conformational resampling from the seeds is partly similar to parallel cascade molecular dynamics (PaCS-MD) simulations, which was developed by one of us.[22] Herein, OFLOOD has two different characteristics against PaCS-MD. (i) OFLOOD only needs reactant(s), although PaCS-MD needs two endpoint structures, that is, a set of a reactant and a product to generate conformational transition pathways connecting these endpoint structures. (ii) In OFLOOD, to enhance the conformational sampling, states 98
Journal of Computational Chemistry 2015, 36, 97–102
Figure 1. Conformational resampling from the outliers (red) promotes the transitions to the neighboring clusters. FlexDice detects outliers by hierarchical clustering in a top-down fashion. In the current case, the twodimensional distributions (black) are divided into three cells at the kth layer, (dense cell: deep green, middle cell: light green, sparse cell: light blue). Then, the divided middle cells at the kth layer are further divided into dense or sparse cells at the k 1 1 layer. In the final stage (the bottom), the divided dense cells are connected between the neighboring dense cells. As shown by arrows, the transitions are stochastic processes depending on the initial velocities in restarting MD simulations. The arrows (blue) represent the cases in which the outliers successfully transited to the neighboring clusters. In contrast, the arrows (purple) represent the cases that the outliers unfortunately failed to transit, but expanded the overall distributions (yellow).
with low population are detected as the outliers by the clustering and the conformational sampling is repeated for the outliers. In this sense, OFLOOD always focuses on the outliers that have low populations in the conformational distributions, some of which have a high potential to across an energy barrier to reach another stable state. In PaCS-MD, structures to be resampled are selected based on several measures with respect to products. Therefore, the structures selected in PaCSMD include states of biomolecules with not only low but also relatively high populations in the conformational distributions.
Results and Discussion As a first demonstration, OFLOOD was applied to extract folding processes of an artificial mini-protein Chignolin (10 residues and 138 atoms),[24] from a completely extended to native structures. To include the solvent effect implicitly, the generalized born with surface area (GBSA) model was used. The salt concentration was set to 0.2 M, and the surface tension was 0.005 kcal/mol/A˚[2] (IGB 5 5 in AMBER 11).[25] MD simulations were performed using the SANDER module in AMBER 11,[25] with the AMBER 99SB force field,[26] which are used in the preceding studies.[27,28] The temperature of the system was controlled at 200 K using the Berendsen thermostat[29] to confirm the conformational sampling efficiency even at a lower temperature. For RCs, the distances of the hydrogen bond (HB1: Asp3N-Gly7O, HB2: Asp3N-Thr8O) were selected, they are essential for the folding of Chignolin.[27,28] The clustering was performed using two-dimensional distribution projected onto WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
Figure 3. Cumulative distributions projected onto each RC with error bars (green) over the last three cycles, a and b) for Chignolin and c and d) for HP35. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
Figure 2. a–d) Projections of the trajectories (red) after the conformational resampling from the outliers (black) at the second, fourth, sixth, and eighth cycles, respectively. The cross in (b) represents the native state of Chignolin. e) The population (the minus common logarithm of frequencies) of the reactive trajectories jointing the snapshots in the initial and final cycles. The minor-folding pathway extracted through OFLOOD is indicated by the dashed arrows. f ) The projections of the trajectories of the 1-ms CMD simulation, initiated from a completely extended structure.
the subspace spanned by HB1 and HB2. To detect the initial outliers, we used a trajectory that was sampled by 10-ns MD simulations started from the completely extended structure. As a reference, a longer (1-ms) CMD simulation was performed from the same initial structure. The detection of outliers and their conformational resampling were continued until 10 cycles. The numbers of detected outliers during the 10 cycles were as follows: (229, 312, 285, 145, 209, 275, 336, 191, 91, and 246). For the conformational resampling, 100-ps MD simulation was restarted from each outlier via regenerations of initial velocities based on the Maxwell Boltzmann distribution. The red points in Figures 2a–2e correspond to the projections of the trajectories after the conformational resampling at the second, fourth, sixth, and eighth cycles, respectively. As shown in Figure 2b, the native state highlighted by the cross, (HB1, ˚ , 3.0 A˚), was sufficiently sampled at the fourth HB2) (6.0 A cycle. By repeating the conformational resampling from the outliers, the projections of the outliers gradually became broader in the subspace. As shown in Figures 2a–2d, the conformational resampling from the outliers filled the unvisited areas of the subspace. In contrast, Figure 2f showed the longer (1-ms) CMD failed to find the native state, leading to the insufficient conformational sampling due to trapping into the mis˚ , 6.0 A˚). fold state near (HB1, HB2) (3.0 A
To monitor the termination of OFLOOD, the convergences of the cumulative distributions projected onto HB1 and HB2 were assessed. Figures 3a–3b show the cumulative distributions averaged over the last three cycles, with the standard deviations. Judging from the error bars averaged over the last three cycles, 10 cycles were sufficient to terminate OFLOOD, as the error bars were quite small. In the next, we discuss the physical meaning of the clusters and the outliers extracted in this study. To classify the stable states, the clustering was performed using all of the trajectories obtained from OFLOOD. As shown in Figure 4a, the stable states appeared as several clusters, and the others were extracted as outliers. Through the outlier analyses, we might obtain rough information about the meta-stable states and their networks interactive. For instance, several folding pathways might be presumed by tracing the distributions of the outliers with the first and second degrees (the black and
Figure 4. a) Projections of the clusters and the outliers obtained from FlexDice. The clusters (red) and the outliers with the first (black), second (yellow), third (cyan), and fourth (magenta) degrees, respectively. (b) The folding FEL at 300 K, obtained from the 2-ms REMD. The values of FEL are scaled by kBT. The major-folding pathway is indicated by the dashed arrows. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
Journal of Computational Chemistry 2015, 36, 97–102
99
FULL PAPER
WWW.C-CHEM.ORG
yellow distributions in Fig. 4a). Each cluster is surrounded by these outliers, indicating that it might exist as transitional states among the clusters. To validate the presumed metastable states and folding pathways, a long-time replica exchange molecular dynamics (REMD) simulation was performed in implicit water. In REMD, 20 replicas (different temperatures distributed from 300 to 490 K, at 10 K intervals) were prepared. The REMD simulation was performed, by exchanging the replicas with adjacent temperatures, until 100 ns were achieved for each replica (2-ls in total). Figure 4b shows the FEL at 300 K obtained from the 2-ms REMD. As Figures 4a–4b show, the distributions of the clusters agreed well with low free energy areas in Figure 4b. The distributions of the outliers with the first and second degrees showed relatively good correspondences with the transitional areas in the FEL of the long-time REMD. Judging from this evidence, the conformational resampling from the outliers tends to promote structural transitions among clusters (meta-stable states), indicating that the outliers might be relevant states for transiting to the neighboring clusters. These simple analyses provide coarse-grained descriptions of the structural transitions. To address the extracted folding pathways, we joined the multiple trajectories in each cycle as reactive trajectories. Herein, reactive trajectories were defined as jointed trajectories connecting the snapshots from the first to last cycles. Figure 2e shows the population of the reactive trajectories projected onto the subspace during the 10 cycles. As an analysis of tracing the outliers, a minor-folding pathway (the dashed arrows in Fig. 2e) that differs from the major one might be presumed, thus indicating that OFLOOD extracted the rare event that directly folds into the native state passing through the misfolded state. The representative snapshots along the minor-folding pathway, which connects the misfolded to the native states, are shown in Figure 5. In the misfolded state, the hydrophobic residues (Tyr2 and Trp9) essential for stabilizing the native state are separated without forming p-p or T-shape stacking due to the hydrophobic interaction. To directly transit to the native state, these hydrophobic residues contacted to make a core (transition state 1) from the misfolded state to undo hydrogen bonds and it was once collapsed (transition state 2) to form the hydrogen bonds suitable to the native state, finally reached the correct T-shape stacking as the “edge to face” hydrophobic interaction. These minor-folding pathways are necessary to obtain a correct statistical ensemble as energetically possible folding processes, despite the fact that they might be unfavorable due to their high free energy profiles. As shown in Figure 4b, the minor-folding pathway was not observed in the FEL at room temperature (300 K) obtained from the 2-ms REMD. In contrast, OFLOOD broadly swept the subspace during the 10 cycles with ns-order computational cost (231.9 ns) even at low temperature (200 K), confirming the high conformational sampling efficiency of OFLOOD, which is a severe test in comparison with a room temperature at 300 K. As a more realistic demonstration, Villin headpiece subdomain (HP35) was considered and its folding processes were investigated by OFLOOD. HP35 is a 35-residue protein consisting of three bundle helices, amounting to 582 atoms. As an initial structure, a completely extended chain was modeled 100
Journal of Computational Chemistry 2015, 36, 97–102
Figure 5. Representative snapshots of the minor-folding pathway connecting the misfolded and the native structures of Chignolin. The characteristic residues for stabilizing the native state, Tyr2 (green) and Trp9 (white), are represented by licorice. The backbone is also represented by ribbon colored with red. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
based on the amino-acid sequence. To compare the conformational sampling efficiency, the same parameters with the previous study (8-ms REMD simulations of HP35) were used.[30] In this demonstration, a set of partial RMSDs that are defined by two segments,[30] segment A (residues 3–21) and segment B (residues 15–33), were selected as RCs. To obtain the outlier at the first cycle, a hundred 100-ps MD simulations with different initial velocities were started from the completely extended structure under NVT (T 5 300 K). Then, FlexDice detected the first outliers and OFLOOD was continued until 20 cycles. We have also confirmed that the cumulative distributions converge within the 20 cycles (see Figs. 3c–3d). The numbers of detected outlier during the 20 cycles were as followings: (99, 96, 82, 126, 73, 75, 66, 77, 101, 126, 79, 76, 88, 74, 66, 52, 40, 53, 61, and 65). The conformational resampling was performed by repeating a set of 100-ps MD simulations started from outliers detected in each cycle. Figures 6a–6e show the projections of trajectories onto the subspace spanned by the partial RMSDs every five cycles, respectively. In these figures, the outliers (black) gradually expanded to low RMSD regions according to the conformational resampling from the outliers. Figure 6 shows the projections of the conformational resampling from the outliers onto the subspace spanned by the partial RMSDs according to the cycles. Here, OFLOOD sufficiently sampled the native structure at the 15th cycle (totally 135.6 ˚ . In Figure 6e, ns) by the following criterion: Ca RMSDs < 1.0 A the structures determined by X-ray (red) and sampled by OFLOOD (blue) are superimposed, showing a good agreement. ˚ . For the folding procHere, the minimum Ca RMSD was 0.5 A esses, OFLOOD might extract the minor-folding pathway in addition to the major-folding pathway. Figure 6f shows the projections of the reactive trajectories during the 20 cycles, in which the reactive trajectories were branched to two folding pathways. The minor-folding pathway has been not sampled in the previous study through the 8-ms REMD simulations,[30] although OFLOOD sampled the rare events with a hundred nsorder computational cost, which also shows the highly conformational sampling efficiency of OFLOOD. Finally, we here would like to mention that OFLOOD is applicable to high-dimensional distributions of states of biomolecules, although we treated simple folding simulations of the mini-proteins in the reduced low-dimensional spaces. For the clustering of the high-dimensional distributions, FlexDice has a robust algorithm to rapidly perform the clustering to detect the outliers. If the high-dimensional distributions can be reduced to WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
Figure 6. Conformational resampling from the outliers (black) at the first, fifth, tenth, fifteenth, and twentieth cycles, respectively. a–e) The crosses in (e and f ) represent the native structure of HP35. f ) Projections of the reactive trajectories of HP35. The major and minor-folding pathways are depicted by the dashed arrows. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
low-dimensional ones, the clustering will be effectively performed in the reduced low-dimensional space spanned by a set of adequate RCs. Without the reductions to low-dimensional distributions, the effectiveness of the clustering seems to be low compared to that of the reduced ones, since the highdimensional distributions tend to be sparse with short-time simulations as performed in this article and difficult to detect the important outliers from them, so that it is necessary to perform relatively long-time simulation in the resampling. In such cases, the most of the components will be regarded as noises. Therefore, it is desirable to reduce the high-dimensional distributions to lower ones using a set of characteristic RCs. In this sense, a choice of RCs directly affects the enhancement of the conformational sampling of OFLOOD, which means that how to specify appropriate RCs is quite important. Actually, the choices of RCs are nontrivial and depend on target systems. Especially, this is a remaining issue in the protein folding simulation and specifications of a set of appropriate RCs to characterize the folding events are still difficult. That is why we simply used a few RCs specified in the past studies. However, the RCs specified in the past studies contain the information of the native structures such like the hydrogen bonds of the backbone (Chignolin) or the partial RMSDs measured from the native structures (HP35),[31] meaning that we implicitly used the information of the native structures in predicting. That is why our predictions of the native structures might be not rigorously blind ones.
FULL PAPER
Therefore, it is instructive to confirm different kind of RCs for the protein folding problems. As a more rigorous and blind prediction, we here would like to show an example of the folding of the same system (Chignolin) without using the information of the native structure. In the new prediction, we used Ca RMSD measured from the initial (unfolded) structure and a radius gyration (Rg) of the all atoms. Based on these newly specified RCs, FlexDice detected outliers from the reduced distributions in the new subspace. Then, the conformational resampling from the new outliers was repeated. In the newly specified RCs, the conformational resampling from the outliers was repeated until 20 cycles. For the computational details on the number of the seeds in each cycle, see the Supporting Information. The projections of the conformational resampling from the outliers are shown every five cycles (Supporting Information Figure S1). In Supporting Information Figure S1, the left column corresponds to the projections onto the subspace spanned by the new RCs. To compare with the previous results, the new trajectories were reprojected onto the original RCs spanned by HB1 and HB2. As shown in Supporting Information Figs. S1(a-b), the native state was still not sampled in the fifth cycle. However, it seems that the native state was sufficiently sampled within the 10th cycle (see Supporting Information Figs. S1(c-d)). Actually, totally 25.3 ns computational time was required to sample the native state using with the new RCs. Therefore, it is numerically demonstrated that OFLOOD could sample the native structure without using the information on the native state in this specific case. As other specifications, there are several candidates of adequate RCs for the protein folding. For instance, the end-to-end distance or the contact number between residues might be appropriate in the protein folding, as the folded structures tend to be globular and compact forms with large number of residue–residue contacts in the native states.[32,33] As the above examples shown, one has some clues to specify RCs in most of applications and can select a set of RCs based on our experiences. Although these RCs were found to be empirically useful for describing the protein folding, more suitable RCs might be carefully considered depending on target system. Herein, we would like to mention several specifications of RCs in applications. For collective movements on domain motions of proteins, principal coordinates defined by the principal component analysis (PCA)[34,35] might be appropriate, as the most of the domain motions are related to anisotropic principle modes (PMs) as low frequency vibrational modes. In the example of the open-closed transition on T4 lysozyme (T4L), the first and the second PMs covered about 70% of the overall fluctuations.[23] Therefore, the high-dimensional distributions will drastically reduced to the low-dimensional distributions in the subspace spanned by a few PMs, which helps us to effectively perform the clustering. As another example, allosteric effects on ligand binding might be difficult to treat by only the PCA, as local and global fluctuations are correlated between the solutes and the ligands in their functions, which means that one have to consider DOF of both solutes and ligands, requiring that suitable RCs to characterize the dynamical correlation between them. As the allosteric effect, we Journal of Computational Chemistry 2015, 36, 97–102
101
FULL PAPER
WWW.C-CHEM.ORG
would like to introduce our recent study about the induced-fit mechanism of substrate enzyme of binding structures of the nylon-origomer hydrolase using PaCS-MD.[36] In this application, a partial RMSD defined by the loop region of this biomolecule, which is relevant region to the allosteric effect on the ligand binding, were specified as RCs. Compared with a case of specification of RC defined by the overall RMSD, the partial one could more effectively promote the structural transitions and were suitable RCs for describing the allosteric effects. Therefore, the choice of RCs depends on target systems and phenomena, and they are sensitive to the conformational sampling efficiency of OFLOOD. How to specify a set of suitable RCs is a future perspective in this study.
Conclusions In this study, we proposed an efficient conformational sampling method, OFLOOD. In this method, the outliers that correspond to sparse distributions in the conformational space were extracted and regarded as relevant states for the structural transitions of proteins. Based on the detections of outliers, we demonstrated that biologically rare events tend to be induced from the outliers through their conformational resampling via restarting short-time MD simulations. In this process, the reassignments of initial velocities might provide sufficient kinetic energies to systems to cross over free energy barriers among meta-stable states. This concept indicates that tracing the outliers leads to extractions of biologically rare events via expansions of highdimensional distributions in the conformational space. From the point of view of applications, OFLOOD might be easy to implement and apply to broad biological targets, once highdimensional distributions projected onto RCs were obtained.
[3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
[15]
[16] [17] [18] [19] [20] [21] [22] [23] [24] [25]
Acknowledgments [26]
R.H thanks Dr. Y. Matsunaga and Dr. Y. Yamamori for helpful discussions. Y.S is thankful to Dr. R. Kishi and Dr. F. Pietrucci for their comments on this work. The calculations were performed using the RIKEN Integrated Cluster of Clusters (RICC) facility. Keywords: biologically rare events structural transitions of proteins rare event sampling outlier detections conformational resampling
How to cite this article: R. Harada, T. Nakmura, Y. Takano, Y. Shigeta. J. Comput. Chem. 2015, 36, 97–102. DOI: 10.1002/ jcc.23773
]
Additional Supporting Information may be found in the online version of this article.
[1] K. Lindorff-Larse, S. Piana, R. O. Dror, D. E. Shaw, Science 2011, 334, 517. [2] R. H. Swendsen, J. S. Wang, Phys. Rev. Lett. 1986, 57, 2607.
102
Journal of Computational Chemistry 2015, 36, 97–102
[27] [28] [29] [30] [31] [32] [33] [34] [35] [36]
K. Hukushima, K. Nemoto, J. Phys. Soc. Jpn. 1996, 65, 1604. Y. Sugita, Y. Okamoto, Chem. Phys. Lett. 1999, 314, 141. B. A. Berg, T. Celik, Int. J. Mod. Phys. C 1992, 03, 1251. A. Mitsutake, Y. Okamoto, J. Chem. Phys. 2009, 130, 214105. A. Mitsutake, J. Chem. Phys. 2009, 131, 094105. B. A. Berg, T. Neuhaus, Phys. Rev. Lett. 1992, 88, 9. N. Nakajima, H. Nakamura, A. Kidera, J. Phys. Chem. B 1997, 101, 817. U. H. E. Hansmann, Y. Okamoto, F. Eisenmenger, Chem. Phys. Lett. 1996, 259, 321. A. Laio, M. Parrinello, Proc. Nat. Acad. Sci. USA 2002, 99, 12562. C. Dellago, P. G. Bolhuis, F. S. Csajka, D. Chandler, J. Chem. Phys. 1998, 108, 1964. W. E, Weiqing Ren, E. Vanden-Eijnden, Chem. Phys. Lett. 2005, 413, 242. T. Nakamura, Y. Kamidoi, S. Wakabayashi, N. Yoshida, Proceedings of the international Special Workshop on Databases for Next Generation Researchers (SWOD2005), Tokyo Japan, April 5-8, 2005, Ed.; IEEE Computer Press. T. Nakamura, Y. Kamidoi, S. Wakabayashi, N. Yoshida, Proceedings of the Research Issues in Data Engineering Stream Data Mining and Applications (RIDE-SDMA 2005), Tokyo Japan, April 3-4, 2005, Ed.; IEE Computer Press. T. Nakamura, Y. Kamidoi, S. Wakabayashi, N. Yoshida, IPSJ Journal Database 2005, 46, 40. T. Nakamura, Y. Kamidoi, S. Wakabayashi, N. Yoshida, IEEE Computer Press. SWOD2006, 2006. L. S. D. Caves, J. D. Evanseck, M. Karplus, Protein Sci. 1998, 7, 649. G. A. Huber, S. Kim, Biophys. J. 1996, 70, 97. B. W. Zhang, D. Jasnow, D. M. Zuckerman, J. Chem. Phys. 2010, 135, 054107. M. Nakano, H. Watanabe, S. M. Rothstein, S. Tanaka, J. Phys. Chem. B 2010, 114, 10234. R. Harada, A. Kitao, J. Chem. Phys. 2013, 139, 035103. R. Harada, Y. Takano, Y. Shigeta, J. Chem. Phys. 2014, 140, 125103. S. Honda, K. Yamasaki, Y. Sawada, H. Morii, Structure 2004, 12, 1507. D. A. Case, T. A. Darden, T. E. Cheatham, III, C. L. Simmerling, J. Wang, R. E. Duke, R., R. C. W. Luo, W. Zhang, K. M. Merz, B. Roberts, B. Wang, S. Hayik, A. Roitberg, I. K. G. Seabra, K. F. Wong, F. Paesani, J. Vanicek, J. Liu, X. Wu, S. R. Brozell, H. G. T. Steinbrecher, Q. Cai, X. Ye, J. Wang, M.- J. Hsieh, G. Cui, D. R. Roe, D. H., M. G. S. Mathews, C. Sagui, V. Babin, T. Luchko, S. Gusarov, A. Kovalenko, P. A. Kollman, AMBER11; University of Carlifornia: San Francisco, 2010. V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg, C. Simmerling, Proteins 2006, 65, 712. D. Satoh, K. Shimizu, S. Nakamura, T. Terada, FEBS. Lett. 2006, 580, 3422. R. Harada, A. Kitao, J. Phys. Chem. B 2011, 115, 8806. H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola, J. R. Haak, J. Chem. Phys. 1984, 81, 3684. L. Hongxing, C. Wu, L. Haiguang, D. Yong, Proc. Nat. Acad. Sci. USA 2007, 104, 4925. T. K. Chiu, J. Kubelka, R. Herbst-Irmer, W. A. Eaton, J. Hofrichter, D. R. Davies, Proc. Natl. Acad. Sci. USA 2005, 102, 7517. F. B. Sheinerman, C. L. Brooks, J. Mol. Biol. 1998, 278, 439. M. Kouza, U. H. E. Hansmann, J. Phys. Chem. B 2012, 116, 6645. A. Amadei, A. B. Linssen, H. J. Berendsen, Proteins 1993, 17, 412. A. Kitao, F. Hirata, N. Go, Chem. Phys. 1991, 158, 447. T. Baba, R. Harada, M. Nakano, Y. Shigeta, J. Comput. Chem. 2014, 35, 1240.
Received: 22 August 2014 Revised: 8 October 2014 Accepted: 16 October 2014 Published online on 3 November 2014
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
Protein Folding Pathways Extracted by OFLOOD: Outlier FLOODing Method Ryuhei Harada,*[a,b] Tomotake Nakamura,[c] Yu Takano,[b,d] and Yasuteru Shigeta[a,b] The Outlier FLOODing method (OFLOOD) is proposed as an efficient conformational sampling method to extract biologically rare events such as protein folding. In OFLOOD, sparse distributions (outliers in the conformational space) were regarded as relevant states for the transitions. Then, the transitions were enhanced through conformational resampling from the outliers. This evidence indicates that the conformational resampling of the sparse distributions might increase chances for promoting the transitions from the outliers to other meta-stable states, which resembles a conformational flooding from the outliers to the neighboring clusters.
OFLOOD consists of (i) detections of outliers from conformational distributions and (ii) conformational resampling from the outliers by molecular dynamics (MD) simulations. Cycles of (i) and (ii) are simply repeated. As demonstrations, OFLOOD was applied to folding of Chignolin and HP35. In both cases, OFLOOD automatically extracted folding pathways from unfolded structures with ns-order computational costs, although ms-order canonical MD failed to extract them. C 2014 Wiley Periodicals, Inc. V
Introduction
used for many targets, such as spin glasses[5] and folding of small peptides.[6,7] Those with the non-Boltzmann weight, such as Multicanonical Monte Carlo[8] and Multicanonical MD simulations,[9,10] were also developed and are commonly used in biophysics for the enhanced conformational sampling of biomolecules. As another powerful structural sampling method, meta-dynamics[11] provides FELs projected onto a subspace spanned by a set of reaction coordinates (RCs) using historydependent repulsive biasing potentials. As strategies for searching the transition pathways, transition path sampling[12] and the string method[13] have rigorous formulations to obtain conformational transition path ensembles. The above conformational sampling methods have been well established and widely used in several biological systems. However, there might exist an issue to be reconsidered when applying these methods to large systems. Generally, with the increase in the vast number of degrees of freedom (DOF), the conformational spaces to be swept become complicated due to high dimensionalities.
Rare events in biological systems are strongly related to their functions, such as ligand binding or allosteric regulation, and so on. In most biological processes, the rare events are observed as structural transitions on biological functions, and thus their extractions are indispensable for understanding biological relevancy. Molecular dynamics (MD) simulation is one of the primary tools to assess biologically important reactions, as a time-series of trajectories at atomic resolution. However, the time scale of biological functions might exceed the accessible time scale of conventional molecular dynamics (CMD) simulations. Recently, a massive supercomputer designed by D. E. Shaw Group has simulated the microsecond-order folding simulations of 12 structurally diverse proteins and unveiled a series of folding mechanisms at atomic levels.[1] From the point of view of computational costs, these long-time CMD simulations seem expensive. In addition, the biologically rare events such as protein folding tend to be observed as stochastic processes. Therefore, it does not ensure that the long-time CMD simulations surely extract the rare events due to the uncertainties coming from the stochastic processes. From the point of the view of free energy landscape (FEL) calculations, CMD simulations with the Boltzmann weight might fail to sample the broad conformational spaces because of trapping into local minima separated by high free energy barriers. To tackle the insufficiency of the structural sampling, developments of methods for promoting the structural transitions are also indispensable for extracting the biologically rare events. Herein, there are two fundamental issues to solve in this field. One is how to find the global and local minima, which are biologically relevant. The other is how to draw FELs. As methods for investigating broad FELs, generalized ensemble algorithms with the Boltzmann weight, such as the parallel tempering[2] and replica exchange method (REM),[3,4] have been widely
DOI: 10.1002/jcc.23773
[a] R. Harada, Y. Shigeta Division of Life Science, Center for Computational Sciences, University of Tsukuba, Tennodai, Tsukuba, Ibaraki 305-8577, Japan E-mail: [email protected] or [email protected] [b] R. Harada, Y. Takano, Y. Shigeta JST-CREST, Kawaguchi, Saitama 332-0012, Japan [c] T. Nakamura Fujitsu Limited, Numazu, Shizuoka 411-0301, Japan [d] Y. Takano Research Center for State-of-the-Art Functional Protein Analysis, Institute for Protein Research, Osaka University, Suita, Osaka 565-0871, Japan Contract grant sponsor: JSPS KAKENHI (Young Scientists (B)); Contract grant number: 26840056; Contract grant sponsor: Scientific Research on Innovative Areas; Contract grant sponsor: MEXT; Contract grant numbers: 26107004 and 26105012; Contract grant sponsor: Computational Materials Science Initiative (CMSI; Consortium) C 2014 Wiley Periodicals, Inc. V
Journal of Computational Chemistry 2015, 36, 97–102
97
FULL PAPER
WWW.C-CHEM.ORG
In treating complicated biological systems, the key point is how to treat the high dimensionalities of the conformational spaces due to increasing of DOF. To describe states of biological systems, a set of RCs should be selected and highdimensional spaces are defined by RCs. Here, to analyze the high-dimensional distributions, clustering techniques might be powerful, as they directly and systematically treat the high dimensionalities. In the clustering processes, stable states of a given protein are regarded as dense distributions in the conformational spaces, that is, clusters. In contrast, sparse distributions are detected as outliers that do not belong to the clusters. These outliers might have potentials to transit to neighboring clusters with high probabilities. Therefore, in this study, we regard the outliers as relevant states that existing in transitional regions. Furthermore, we show that the conformational resampling from the outliers via short-time MD simulations tends to promote the structural transitions of proteins, which efficiently leads to extractions of biologically rare events.
Methods Based on the detections of outliers and their conformational resampling, we propose an efficient conformational sampling method to extract the biologically rare events of proteins, called Outlier FLOODing method (OFLOOD). In OFLOOD, a set of cycle consisting of the following schemes is simply repeated: (i) detections of the outliers from conformational distributions by clustering of MD trajectories and (ii) conformational resampling from the detected outliers through shorttime MD simulations. Here, the detections of the outliers and clustering of MD trajectories are performed using a newly developed clustering algorithm, FlexDice.[14–17] The conformational resampling is performed by restarting short-time MD simulations from the outliers via regenerations of initial velocities. In several past studies,[18–23] restarting MD simulations via regenerations of initial velocities might enhance the conformational sampling of proteins. In every cycle, outliers are updated based on cumulative MD trajectories obtained from the past conformational resampling. These cycles are continued until the conformational search sufficiently converges, which means that the convergences of cumulative distributions projected onto each RC are checked for the termination of OFLOOD. Figure 1 depicts the physical meaning of tracing the outliers in OFLOOD. As shown in Figure 1, the conformational resampling from the outliers might expand distributions through transitions to the neighboring clusters, which leads to extractions of biologically rare events of proteins. The strategy of OFLOOD based on the conformational resampling from the seeds is partly similar to parallel cascade molecular dynamics (PaCS-MD) simulations, which was developed by one of us.[22] Herein, OFLOOD has two different characteristics against PaCS-MD. (i) OFLOOD only needs reactant(s), although PaCS-MD needs two endpoint structures, that is, a set of a reactant and a product to generate conformational transition pathways connecting these endpoint structures. (ii) In OFLOOD, to enhance the conformational sampling, states 98
Journal of Computational Chemistry 2015, 36, 97–102
Figure 1. Conformational resampling from the outliers (red) promotes the transitions to the neighboring clusters. FlexDice detects outliers by hierarchical clustering in a top-down fashion. In the current case, the twodimensional distributions (black) are divided into three cells at the kth layer, (dense cell: deep green, middle cell: light green, sparse cell: light blue). Then, the divided middle cells at the kth layer are further divided into dense or sparse cells at the k 1 1 layer. In the final stage (the bottom), the divided dense cells are connected between the neighboring dense cells. As shown by arrows, the transitions are stochastic processes depending on the initial velocities in restarting MD simulations. The arrows (blue) represent the cases in which the outliers successfully transited to the neighboring clusters. In contrast, the arrows (purple) represent the cases that the outliers unfortunately failed to transit, but expanded the overall distributions (yellow).
with low population are detected as the outliers by the clustering and the conformational sampling is repeated for the outliers. In this sense, OFLOOD always focuses on the outliers that have low populations in the conformational distributions, some of which have a high potential to across an energy barrier to reach another stable state. In PaCS-MD, structures to be resampled are selected based on several measures with respect to products. Therefore, the structures selected in PaCSMD include states of biomolecules with not only low but also relatively high populations in the conformational distributions.
Results and Discussion As a first demonstration, OFLOOD was applied to extract folding processes of an artificial mini-protein Chignolin (10 residues and 138 atoms),[24] from a completely extended to native structures. To include the solvent effect implicitly, the generalized born with surface area (GBSA) model was used. The salt concentration was set to 0.2 M, and the surface tension was 0.005 kcal/mol/A˚[2] (IGB 5 5 in AMBER 11).[25] MD simulations were performed using the SANDER module in AMBER 11,[25] with the AMBER 99SB force field,[26] which are used in the preceding studies.[27,28] The temperature of the system was controlled at 200 K using the Berendsen thermostat[29] to confirm the conformational sampling efficiency even at a lower temperature. For RCs, the distances of the hydrogen bond (HB1: Asp3N-Gly7O, HB2: Asp3N-Thr8O) were selected, they are essential for the folding of Chignolin.[27,28] The clustering was performed using two-dimensional distribution projected onto WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
Figure 3. Cumulative distributions projected onto each RC with error bars (green) over the last three cycles, a and b) for Chignolin and c and d) for HP35. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
Figure 2. a–d) Projections of the trajectories (red) after the conformational resampling from the outliers (black) at the second, fourth, sixth, and eighth cycles, respectively. The cross in (b) represents the native state of Chignolin. e) The population (the minus common logarithm of frequencies) of the reactive trajectories jointing the snapshots in the initial and final cycles. The minor-folding pathway extracted through OFLOOD is indicated by the dashed arrows. f ) The projections of the trajectories of the 1-ms CMD simulation, initiated from a completely extended structure.
the subspace spanned by HB1 and HB2. To detect the initial outliers, we used a trajectory that was sampled by 10-ns MD simulations started from the completely extended structure. As a reference, a longer (1-ms) CMD simulation was performed from the same initial structure. The detection of outliers and their conformational resampling were continued until 10 cycles. The numbers of detected outliers during the 10 cycles were as follows: (229, 312, 285, 145, 209, 275, 336, 191, 91, and 246). For the conformational resampling, 100-ps MD simulation was restarted from each outlier via regenerations of initial velocities based on the Maxwell Boltzmann distribution. The red points in Figures 2a–2e correspond to the projections of the trajectories after the conformational resampling at the second, fourth, sixth, and eighth cycles, respectively. As shown in Figure 2b, the native state highlighted by the cross, (HB1, ˚ , 3.0 A˚), was sufficiently sampled at the fourth HB2) (6.0 A cycle. By repeating the conformational resampling from the outliers, the projections of the outliers gradually became broader in the subspace. As shown in Figures 2a–2d, the conformational resampling from the outliers filled the unvisited areas of the subspace. In contrast, Figure 2f showed the longer (1-ms) CMD failed to find the native state, leading to the insufficient conformational sampling due to trapping into the mis˚ , 6.0 A˚). fold state near (HB1, HB2) (3.0 A
To monitor the termination of OFLOOD, the convergences of the cumulative distributions projected onto HB1 and HB2 were assessed. Figures 3a–3b show the cumulative distributions averaged over the last three cycles, with the standard deviations. Judging from the error bars averaged over the last three cycles, 10 cycles were sufficient to terminate OFLOOD, as the error bars were quite small. In the next, we discuss the physical meaning of the clusters and the outliers extracted in this study. To classify the stable states, the clustering was performed using all of the trajectories obtained from OFLOOD. As shown in Figure 4a, the stable states appeared as several clusters, and the others were extracted as outliers. Through the outlier analyses, we might obtain rough information about the meta-stable states and their networks interactive. For instance, several folding pathways might be presumed by tracing the distributions of the outliers with the first and second degrees (the black and
Figure 4. a) Projections of the clusters and the outliers obtained from FlexDice. The clusters (red) and the outliers with the first (black), second (yellow), third (cyan), and fourth (magenta) degrees, respectively. (b) The folding FEL at 300 K, obtained from the 2-ms REMD. The values of FEL are scaled by kBT. The major-folding pathway is indicated by the dashed arrows. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
Journal of Computational Chemistry 2015, 36, 97–102
99
FULL PAPER
WWW.C-CHEM.ORG
yellow distributions in Fig. 4a). Each cluster is surrounded by these outliers, indicating that it might exist as transitional states among the clusters. To validate the presumed metastable states and folding pathways, a long-time replica exchange molecular dynamics (REMD) simulation was performed in implicit water. In REMD, 20 replicas (different temperatures distributed from 300 to 490 K, at 10 K intervals) were prepared. The REMD simulation was performed, by exchanging the replicas with adjacent temperatures, until 100 ns were achieved for each replica (2-ls in total). Figure 4b shows the FEL at 300 K obtained from the 2-ms REMD. As Figures 4a–4b show, the distributions of the clusters agreed well with low free energy areas in Figure 4b. The distributions of the outliers with the first and second degrees showed relatively good correspondences with the transitional areas in the FEL of the long-time REMD. Judging from this evidence, the conformational resampling from the outliers tends to promote structural transitions among clusters (meta-stable states), indicating that the outliers might be relevant states for transiting to the neighboring clusters. These simple analyses provide coarse-grained descriptions of the structural transitions. To address the extracted folding pathways, we joined the multiple trajectories in each cycle as reactive trajectories. Herein, reactive trajectories were defined as jointed trajectories connecting the snapshots from the first to last cycles. Figure 2e shows the population of the reactive trajectories projected onto the subspace during the 10 cycles. As an analysis of tracing the outliers, a minor-folding pathway (the dashed arrows in Fig. 2e) that differs from the major one might be presumed, thus indicating that OFLOOD extracted the rare event that directly folds into the native state passing through the misfolded state. The representative snapshots along the minor-folding pathway, which connects the misfolded to the native states, are shown in Figure 5. In the misfolded state, the hydrophobic residues (Tyr2 and Trp9) essential for stabilizing the native state are separated without forming p-p or T-shape stacking due to the hydrophobic interaction. To directly transit to the native state, these hydrophobic residues contacted to make a core (transition state 1) from the misfolded state to undo hydrogen bonds and it was once collapsed (transition state 2) to form the hydrogen bonds suitable to the native state, finally reached the correct T-shape stacking as the “edge to face” hydrophobic interaction. These minor-folding pathways are necessary to obtain a correct statistical ensemble as energetically possible folding processes, despite the fact that they might be unfavorable due to their high free energy profiles. As shown in Figure 4b, the minor-folding pathway was not observed in the FEL at room temperature (300 K) obtained from the 2-ms REMD. In contrast, OFLOOD broadly swept the subspace during the 10 cycles with ns-order computational cost (231.9 ns) even at low temperature (200 K), confirming the high conformational sampling efficiency of OFLOOD, which is a severe test in comparison with a room temperature at 300 K. As a more realistic demonstration, Villin headpiece subdomain (HP35) was considered and its folding processes were investigated by OFLOOD. HP35 is a 35-residue protein consisting of three bundle helices, amounting to 582 atoms. As an initial structure, a completely extended chain was modeled 100
Journal of Computational Chemistry 2015, 36, 97–102
Figure 5. Representative snapshots of the minor-folding pathway connecting the misfolded and the native structures of Chignolin. The characteristic residues for stabilizing the native state, Tyr2 (green) and Trp9 (white), are represented by licorice. The backbone is also represented by ribbon colored with red. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
based on the amino-acid sequence. To compare the conformational sampling efficiency, the same parameters with the previous study (8-ms REMD simulations of HP35) were used.[30] In this demonstration, a set of partial RMSDs that are defined by two segments,[30] segment A (residues 3–21) and segment B (residues 15–33), were selected as RCs. To obtain the outlier at the first cycle, a hundred 100-ps MD simulations with different initial velocities were started from the completely extended structure under NVT (T 5 300 K). Then, FlexDice detected the first outliers and OFLOOD was continued until 20 cycles. We have also confirmed that the cumulative distributions converge within the 20 cycles (see Figs. 3c–3d). The numbers of detected outlier during the 20 cycles were as followings: (99, 96, 82, 126, 73, 75, 66, 77, 101, 126, 79, 76, 88, 74, 66, 52, 40, 53, 61, and 65). The conformational resampling was performed by repeating a set of 100-ps MD simulations started from outliers detected in each cycle. Figures 6a–6e show the projections of trajectories onto the subspace spanned by the partial RMSDs every five cycles, respectively. In these figures, the outliers (black) gradually expanded to low RMSD regions according to the conformational resampling from the outliers. Figure 6 shows the projections of the conformational resampling from the outliers onto the subspace spanned by the partial RMSDs according to the cycles. Here, OFLOOD sufficiently sampled the native structure at the 15th cycle (totally 135.6 ˚ . In Figure 6e, ns) by the following criterion: Ca RMSDs < 1.0 A the structures determined by X-ray (red) and sampled by OFLOOD (blue) are superimposed, showing a good agreement. ˚ . For the folding procHere, the minimum Ca RMSD was 0.5 A esses, OFLOOD might extract the minor-folding pathway in addition to the major-folding pathway. Figure 6f shows the projections of the reactive trajectories during the 20 cycles, in which the reactive trajectories were branched to two folding pathways. The minor-folding pathway has been not sampled in the previous study through the 8-ms REMD simulations,[30] although OFLOOD sampled the rare events with a hundred nsorder computational cost, which also shows the highly conformational sampling efficiency of OFLOOD. Finally, we here would like to mention that OFLOOD is applicable to high-dimensional distributions of states of biomolecules, although we treated simple folding simulations of the mini-proteins in the reduced low-dimensional spaces. For the clustering of the high-dimensional distributions, FlexDice has a robust algorithm to rapidly perform the clustering to detect the outliers. If the high-dimensional distributions can be reduced to WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
Figure 6. Conformational resampling from the outliers (black) at the first, fifth, tenth, fifteenth, and twentieth cycles, respectively. a–e) The crosses in (e and f ) represent the native structure of HP35. f ) Projections of the reactive trajectories of HP35. The major and minor-folding pathways are depicted by the dashed arrows. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
low-dimensional ones, the clustering will be effectively performed in the reduced low-dimensional space spanned by a set of adequate RCs. Without the reductions to low-dimensional distributions, the effectiveness of the clustering seems to be low compared to that of the reduced ones, since the highdimensional distributions tend to be sparse with short-time simulations as performed in this article and difficult to detect the important outliers from them, so that it is necessary to perform relatively long-time simulation in the resampling. In such cases, the most of the components will be regarded as noises. Therefore, it is desirable to reduce the high-dimensional distributions to lower ones using a set of characteristic RCs. In this sense, a choice of RCs directly affects the enhancement of the conformational sampling of OFLOOD, which means that how to specify appropriate RCs is quite important. Actually, the choices of RCs are nontrivial and depend on target systems. Especially, this is a remaining issue in the protein folding simulation and specifications of a set of appropriate RCs to characterize the folding events are still difficult. That is why we simply used a few RCs specified in the past studies. However, the RCs specified in the past studies contain the information of the native structures such like the hydrogen bonds of the backbone (Chignolin) or the partial RMSDs measured from the native structures (HP35),[31] meaning that we implicitly used the information of the native structures in predicting. That is why our predictions of the native structures might be not rigorously blind ones.
FULL PAPER
Therefore, it is instructive to confirm different kind of RCs for the protein folding problems. As a more rigorous and blind prediction, we here would like to show an example of the folding of the same system (Chignolin) without using the information of the native structure. In the new prediction, we used Ca RMSD measured from the initial (unfolded) structure and a radius gyration (Rg) of the all atoms. Based on these newly specified RCs, FlexDice detected outliers from the reduced distributions in the new subspace. Then, the conformational resampling from the new outliers was repeated. In the newly specified RCs, the conformational resampling from the outliers was repeated until 20 cycles. For the computational details on the number of the seeds in each cycle, see the Supporting Information. The projections of the conformational resampling from the outliers are shown every five cycles (Supporting Information Figure S1). In Supporting Information Figure S1, the left column corresponds to the projections onto the subspace spanned by the new RCs. To compare with the previous results, the new trajectories were reprojected onto the original RCs spanned by HB1 and HB2. As shown in Supporting Information Figs. S1(a-b), the native state was still not sampled in the fifth cycle. However, it seems that the native state was sufficiently sampled within the 10th cycle (see Supporting Information Figs. S1(c-d)). Actually, totally 25.3 ns computational time was required to sample the native state using with the new RCs. Therefore, it is numerically demonstrated that OFLOOD could sample the native structure without using the information on the native state in this specific case. As other specifications, there are several candidates of adequate RCs for the protein folding. For instance, the end-to-end distance or the contact number between residues might be appropriate in the protein folding, as the folded structures tend to be globular and compact forms with large number of residue–residue contacts in the native states.[32,33] As the above examples shown, one has some clues to specify RCs in most of applications and can select a set of RCs based on our experiences. Although these RCs were found to be empirically useful for describing the protein folding, more suitable RCs might be carefully considered depending on target system. Herein, we would like to mention several specifications of RCs in applications. For collective movements on domain motions of proteins, principal coordinates defined by the principal component analysis (PCA)[34,35] might be appropriate, as the most of the domain motions are related to anisotropic principle modes (PMs) as low frequency vibrational modes. In the example of the open-closed transition on T4 lysozyme (T4L), the first and the second PMs covered about 70% of the overall fluctuations.[23] Therefore, the high-dimensional distributions will drastically reduced to the low-dimensional distributions in the subspace spanned by a few PMs, which helps us to effectively perform the clustering. As another example, allosteric effects on ligand binding might be difficult to treat by only the PCA, as local and global fluctuations are correlated between the solutes and the ligands in their functions, which means that one have to consider DOF of both solutes and ligands, requiring that suitable RCs to characterize the dynamical correlation between them. As the allosteric effect, we Journal of Computational Chemistry 2015, 36, 97–102
101
FULL PAPER
WWW.C-CHEM.ORG
would like to introduce our recent study about the induced-fit mechanism of substrate enzyme of binding structures of the nylon-origomer hydrolase using PaCS-MD.[36] In this application, a partial RMSD defined by the loop region of this biomolecule, which is relevant region to the allosteric effect on the ligand binding, were specified as RCs. Compared with a case of specification of RC defined by the overall RMSD, the partial one could more effectively promote the structural transitions and were suitable RCs for describing the allosteric effects. Therefore, the choice of RCs depends on target systems and phenomena, and they are sensitive to the conformational sampling efficiency of OFLOOD. How to specify a set of suitable RCs is a future perspective in this study.
Conclusions In this study, we proposed an efficient conformational sampling method, OFLOOD. In this method, the outliers that correspond to sparse distributions in the conformational space were extracted and regarded as relevant states for the structural transitions of proteins. Based on the detections of outliers, we demonstrated that biologically rare events tend to be induced from the outliers through their conformational resampling via restarting short-time MD simulations. In this process, the reassignments of initial velocities might provide sufficient kinetic energies to systems to cross over free energy barriers among meta-stable states. This concept indicates that tracing the outliers leads to extractions of biologically rare events via expansions of highdimensional distributions in the conformational space. From the point of view of applications, OFLOOD might be easy to implement and apply to broad biological targets, once highdimensional distributions projected onto RCs were obtained.
[3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
[15]
[16] [17] [18] [19] [20] [21] [22] [23] [24] [25]
Acknowledgments [26]
R.H thanks Dr. Y. Matsunaga and Dr. Y. Yamamori for helpful discussions. Y.S is thankful to Dr. R. Kishi and Dr. F. Pietrucci for their comments on this work. The calculations were performed using the RIKEN Integrated Cluster of Clusters (RICC) facility. Keywords: biologically rare events structural transitions of proteins rare event sampling outlier detections conformational resampling
How to cite this article: R. Harada, T. Nakmura, Y. Takano, Y. Shigeta. J. Comput. Chem. 2015, 36, 97–102. DOI: 10.1002/ jcc.23773
]
Additional Supporting Information may be found in the online version of this article.
[1] K. Lindorff-Larse, S. Piana, R. O. Dror, D. E. Shaw, Science 2011, 334, 517. [2] R. H. Swendsen, J. S. Wang, Phys. Rev. Lett. 1986, 57, 2607.
102
Journal of Computational Chemistry 2015, 36, 97–102
[27] [28] [29] [30] [31] [32] [33] [34] [35] [36]
K. Hukushima, K. Nemoto, J. Phys. Soc. Jpn. 1996, 65, 1604. Y. Sugita, Y. Okamoto, Chem. Phys. Lett. 1999, 314, 141. B. A. Berg, T. Celik, Int. J. Mod. Phys. C 1992, 03, 1251. A. Mitsutake, Y. Okamoto, J. Chem. Phys. 2009, 130, 214105. A. Mitsutake, J. Chem. Phys. 2009, 131, 094105. B. A. Berg, T. Neuhaus, Phys. Rev. Lett. 1992, 88, 9. N. Nakajima, H. Nakamura, A. Kidera, J. Phys. Chem. B 1997, 101, 817. U. H. E. Hansmann, Y. Okamoto, F. Eisenmenger, Chem. Phys. Lett. 1996, 259, 321. A. Laio, M. Parrinello, Proc. Nat. Acad. Sci. USA 2002, 99, 12562. C. Dellago, P. G. Bolhuis, F. S. Csajka, D. Chandler, J. Chem. Phys. 1998, 108, 1964. W. E, Weiqing Ren, E. Vanden-Eijnden, Chem. Phys. Lett. 2005, 413, 242. T. Nakamura, Y. Kamidoi, S. Wakabayashi, N. Yoshida, Proceedings of the international Special Workshop on Databases for Next Generation Researchers (SWOD2005), Tokyo Japan, April 5-8, 2005, Ed.; IEEE Computer Press. T. Nakamura, Y. Kamidoi, S. Wakabayashi, N. Yoshida, Proceedings of the Research Issues in Data Engineering Stream Data Mining and Applications (RIDE-SDMA 2005), Tokyo Japan, April 3-4, 2005, Ed.; IEE Computer Press. T. Nakamura, Y. Kamidoi, S. Wakabayashi, N. Yoshida, IPSJ Journal Database 2005, 46, 40. T. Nakamura, Y. Kamidoi, S. Wakabayashi, N. Yoshida, IEEE Computer Press. SWOD2006, 2006. L. S. D. Caves, J. D. Evanseck, M. Karplus, Protein Sci. 1998, 7, 649. G. A. Huber, S. Kim, Biophys. J. 1996, 70, 97. B. W. Zhang, D. Jasnow, D. M. Zuckerman, J. Chem. Phys. 2010, 135, 054107. M. Nakano, H. Watanabe, S. M. Rothstein, S. Tanaka, J. Phys. Chem. B 2010, 114, 10234. R. Harada, A. Kitao, J. Chem. Phys. 2013, 139, 035103. R. Harada, Y. Takano, Y. Shigeta, J. Chem. Phys. 2014, 140, 125103. S. Honda, K. Yamasaki, Y. Sawada, H. Morii, Structure 2004, 12, 1507. D. A. Case, T. A. Darden, T. E. Cheatham, III, C. L. Simmerling, J. Wang, R. E. Duke, R., R. C. W. Luo, W. Zhang, K. M. Merz, B. Roberts, B. Wang, S. Hayik, A. Roitberg, I. K. G. Seabra, K. F. Wong, F. Paesani, J. Vanicek, J. Liu, X. Wu, S. R. Brozell, H. G. T. Steinbrecher, Q. Cai, X. Ye, J. Wang, M.- J. Hsieh, G. Cui, D. R. Roe, D. H., M. G. S. Mathews, C. Sagui, V. Babin, T. Luchko, S. Gusarov, A. Kovalenko, P. A. Kollman, AMBER11; University of Carlifornia: San Francisco, 2010. V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg, C. Simmerling, Proteins 2006, 65, 712. D. Satoh, K. Shimizu, S. Nakamura, T. Terada, FEBS. Lett. 2006, 580, 3422. R. Harada, A. Kitao, J. Phys. Chem. B 2011, 115, 8806. H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola, J. R. Haak, J. Chem. Phys. 1984, 81, 3684. L. Hongxing, C. Wu, L. Haiguang, D. Yong, Proc. Nat. Acad. Sci. USA 2007, 104, 4925. T. K. Chiu, J. Kubelka, R. Herbst-Irmer, W. A. Eaton, J. Hofrichter, D. R. Davies, Proc. Natl. Acad. Sci. USA 2005, 102, 7517. F. B. Sheinerman, C. L. Brooks, J. Mol. Biol. 1998, 278, 439. M. Kouza, U. H. E. Hansmann, J. Phys. Chem. B 2012, 116, 6645. A. Amadei, A. B. Linssen, H. J. Berendsen, Proteins 1993, 17, 412. A. Kitao, F. Hirata, N. Go, Chem. Phys. 1991, 158, 447. T. Baba, R. Harada, M. Nakano, Y. Shigeta, J. Comput. Chem. 2014, 35, 1240.
Received: 22 August 2014 Revised: 8 October 2014 Accepted: 16 October 2014 Published online on 3 November 2014
WWW.CHEMISTRYVIEWS.COM
