Dimensionality of Collective Variables for Describing Conformational

Apr 6, 2016 - (20, 21) For MD simulations and the string method calculations of AdK, a Cα-based coarse-grained (CG) model, which we call the domain m...
1 downloads 11 Views 2MB Size
Download PDF
Letter pubs.acs.org/JPCL

Dimensionality of Collective Variables for Describing Conformational Changes of a Multi-Domain Protein Yasuhiro Matsunaga,† Yasuaki Komuro,‡,⊥ Chigusa Kobayashi,† Jaewoon Jung,†,§ Takaharu Mori,‡,§ and Yuji Sugita*,†,‡,§,∥ †

RIKEN Advanced Institute for Computational Science, 7-1-26 minatojima-minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan RIKEN, Theoretical Molecular Science Laboratory, 2-1 Hirosawa, Wako-shi, Saitama 351-0198, Japan § RIKEN, iTHES, 2-1 Hirosawa, Wako-shi, Saitama 351-0198, Japan ∥ RIKEN Quantitative Biology Center, Integrated Innovation Building 7F, 6-7-1 minatojima-minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan ‡

ABSTRACT: Collective variables (CVs) are often used in molecular dynamics simulations based on enhanced sampling algorithms to investigate large conformational changes of a protein. The choice of CVs in these simulations is essential because it affects simulation results and impacts the freeenergy profile, the minimum free-energy pathway (MFEP), and the transition-state structure. Here we examine how many CVs are required to capture the correct transition-state structure during the open-to-close motion of adenylate kinase using a coarse-grained model in the mean forces string method to search the MFEP. Various numbers of large amplitude principal components are tested as CVs in the simulations. The incorporation of local coordinates into CVs, which is possible in higher dimensional CV spaces, is important for capturing a reliable MFEP. The Bayesian measure proposed by Best and Hummer is sensitive to the choice of CVs, showing sharp peaks when the transition-state structure is captured. We thus evaluate the required number of CVs needed in enhanced sampling simulations for describing protein conformational changes.

L

the molecular mechanism for the conformational transitions.6,12 In protein folding, the fraction of native contacts forms a good choice of CVs.13 Piana et al. showed, from their millisecond allatom simulations, that the fraction of native contacts is sufficient to capture intermediate states during folding even for a slow-folding protein like ubiquitin.14 Best et al. recently analyzed folding simulation trajectories of various proteins15 and concluded that nonnative contacts are unimportant in determining folding mechanism. For large conformational changes of proteins, such as functional open-to-close motions of multidomain proteins, several choices of CVs have been proposed and demonstrated using enhanced sampling methods. For example, Abrams and Vanden-Eijnden16 successfully sampled conformational changes of GroEL subunit and HIV-1 gp120 by using the temperature -accelerated molecular dynamics (TAMD).10 In their simulation, Cartesian coordinates of centers of contiguous subdomains were chosen as for CVs, composed of 9 subdomains for GroEL and 14 subdomains for gp120. Vashisth and Brooks demonstrated that applying a potential energy bias in the direction of displacement known from the crystal

arge-scale conformational changes of proteins are generally important for their function.1 The conformational changes involve barrier-crossing transitions on the complex free-energy surfaces, which are difficult to sample in brute-force all-atom molecular dynamics (MD) simulations. Although specialpurpose computers2−4 have extended the time scale of the MD simulations toward hundreds of microseconds or longer, efficient sampling algorithms are still useful to simulate such rare events, especially on general-purpose computers. In the last 20 years, various types of enhanced sampling methods have been proposed for simulating protein folding and large conformational changes.5−10 Many of them require predefined collective variables (CVs) for introducing biased potentials/ forces or resampling techniques. For example, metadynamics5 adds Gaussian-type potentials in the CV space to the original potential energy. The mean forces string method6,11 searches the minimum free energy path (MFEP) in the CV space by imposing restraint potentials to the CVs. The weighted ensemble method9 resamples snapshots according to probabilities of binned cells in the CV space. In the enhanced simulations, the choice of CVs can influence the conformational sampling efficiency as well as the simulation results. Slow degrees of freedom orthogonal to the CV space can slow down the convergence of sampling. A free-energy profile evaluated with poor CV choice may fail to capture transition-state structures, which leads to a misunderstanding of © XXXX American Chemical Society

Received: February 12, 2016 Accepted: April 4, 2016

1446

DOI: 10.1021/acs.jpclett.6b00317 J. Phys. Chem. Lett. 2016, 7, 1446−1451

Letter

The Journal of Physical Chemistry Letters

Figure 1. (A) Free-energy landscape on the first (horizontal axis) and second (vertical axis) PCs of adenylate kinase. Lines indicate MFEPs calculated in 2D (dark blue), 3D (light blue), 10D (yellow), and 20D (red) PC spaces. (B) Projections of the MFEPs of various dimensional PC spaces projected onto the distances between domains’ centers of mass. Representative structures in open, semiclosed, and closed states are shown with tube representation.

structures facilitates functional motions in TAMD simulations.17 One of the major difficulties in choosing CVs is how to reduce their number from the large degree of freedom of the target protein while ensuring both sampling efficiency and the accuracy of transition pathways in the CV space. In particular, while plausible CV choices for sampling efficiency have been well studied, the impact of CV choice on the accuracy of transition pathways and transition states is only marginally understood. Maragliano et al. showed, for accurate description of conformational transition of alanine-dipeptide, that four dihedral angles (rather than two) are required by evaluating committer probabilities.6 Recently, Pan et al. showed that, through the comparison with brute-force simulations by a special-purpose computer Anton, root-mean-square deviation (RMSD) of the flexible region is a better choice for accurate description of transition pathways of EGFR kinase rather than RMSD of the entire protein.8 Principal component analysis (PCA) has been used to describe conformational fluctuations of native protein conformations.18 Early simulation studies19 showed that a few large amplitude principal components (PCs) are sufficient to capture anharmonic dimensions of total conformational fluctuations; however, it remains elusive as to how many PCs are required to capture a conformational transition pathway and the transitionstate structure with high accuracy. Here we apply the mean forces string method6 to find the MFEP that represents the most probable conformational transition pathway. By changing the number of PCs used to define the CV space in the simulations, we systematically examine the quality of the MFEP for capturing the transition-state structure. We tested the open-to-close transition of a small globular protein, adenylate kinase (AdK)1,20−27 in this study. AdK has 214 residues and consists of three relatively rigid domains: the CORE domain (residues 1−29, 68−117, and 161−214), the NMP domain (residues 30−67), and the LID domain (residues 118−167). AdK catalyzes a transphosphorylation, ATP + AMP ↔ 2ADP, and adopts a closed form with bound ligands and an open form in the absence of ligands. This conformational transition is the rate-limiting step in the enzymatic reaction cycle.20,21 For MD simulations and the string method calculations of AdK, a Cα-based coarse-grained (CG) model, which we call the domain motion enhanced (DoME) model,28

is used. This model, recently developed by the authors, has reasonably low free-energy barriers between the open, intermediate (semiclosed), and closed states.28 The semiclosed form, in which only LID takes the closed conformation, has been observed in previous CG MD simulations.24,25,27 In protein folding studies, the quality of 1D reaction coordinates for describing folding and unfolding events has been successfully evaluated by the Bayesian measure devised by Best and Hummer.13,29 Here, by using brute-force simulation trajectory, we apply the same measure to MFEPs searched in the mean forces string method with multidimensional PC spaces and evaluate the quality of the MFEPs for describing the conformational changes. To this end, we introduce progress coordinates r(x) for pathways developed in the context of metadynamics simulation.30 We first project brute-force simulation trajectory x onto multidimensional PC spaces and define the progress coordinate r(x) using projection z30,31 2

M

r (x ) =

∑m = 1 me−λ || z(x) − z(m) || M

2

∑m = 1 e−λ || z(x) − z(m) ||

(1)

where M is the number of images {z(m)} for discretization of a pathway. λ is comparable to the inverse of mean square displacement between successive images, defined as λ = 2.3/ 2 31 [∑M−1 This choice of λ make m=1 ∥z(m + 1)−z(m)∥/(M − 1)] . 30 r(x) change smoothly from 1 to M. By using r(x), the freeenergy profile along the pathway is defined as F(r) = −kBT log peq(r), where peq(r) is an equilibrium density of r. A transition path (TP) is defined as the portion of the trajectory x that exits from the reactant region A and reaches the product region B without crossing back into A. Following Best and Hummer,13 two probabilities are introduced: p(TP|r) is the probability for being on a TP, given that the system is in r, and p(TP) is the fraction of time spent in TPs. According to Bayes’ rule, p(TP|r) is calculated as p(TP|r ) =

p(r |TP)p(TP) peq (r )

(2)

where p(r|TP) is a probability density of r given that the system is in TPs. Around the transition state in r, the probability of observing TPs should be maximal because the transition state is 1447

DOI: 10.1021/acs.jpclett.6b00317 J. Phys. Chem. Lett. 2016, 7, 1446−1451

Letter

The Journal of Physical Chemistry Letters

Figure 2. (A) Free-energy profiles along the MFEPs calculated in the string method simulations with various dimensional PC spaces. Image index is assigned from the open (first image) to the closed state (64th image). Lines are colored from dark blue to dark brown according to the dimensionality of PC space. (B) Bayesian measure p(TP|r) for open-to-semiclosed transition with various dimensional PC spaces. (C). Bayesian measure p(TP|r) for semiclosed-to-closed transition. (D) Comparison of the structures at the transition state (41th image) constructed from 3D (light colors; pale green, pink, and cyan) and 25D (dark colors; green, red, and blue) PCs. The black circle indicates that the 2D pathway fails to describe the intradomain bending of the NMP domain. The lines between the LID and NMP domains correspond to interactive pairs between two domains in the DoME model.

a dividing surface that separates reactant and product. In particular, a good choice of CVs provides a sharp high peak of p(TP|r), indicating infrequent recrossings around the identified transition state.13,29 Figure 1A shows the free-energy landscape of AdK on the 2D PC space spanned with the first and second (largest amplitude) PCs. The landscape was observed by a long brute-force MD simulation. The MFEPs calculated in the mean forces string method with various PC dimensions are superimposed. Regardless of the PC dimensions, MFEPs capture lower freeenergy regions on the PC space, passing through the open, intermediate, and closed states. In Figure 1B, the MFEPs are projected onto geometric coordinate space, the center-of-mass distances between the CORE and NMP and those between the CORE and LID. Only the MFEP with the 2D PC deviates from other MFEPs, indicating that it fails to capture the probable transition pathway of AdK. The MFEP with the 3D PC is very similar to those with higher dimensional PCs. The first and second PC modes represent correlated motions of the LID and NMP domains, whereas the third PC mode captures the uncorrelated motion of the two domains. Thus, independent motions of the LID and NMP domains, such as the LID fastclosing motion, cannot be represented only by the first and second PC modes. Figure 2A shows the free-energy profile along the MFEPs with various dimensional PC (from the 2D to 25D PCs). Image index is assigned from the open (first image) to the closed state (64th image) by using the progress coordinate r(x) for each MFEP. The semiclosed state is seen in the 31th image. Two

prominent free-energy barriers between the open and semiclosed state and the semiclosed to closed state are observed. Interestingly, as the PC dimension increases the height of the barriers also gradually increases, while the relative stabilities of the three stable states do not change. These increases start to converge when we apply the 10D PC or higher spaces in the mean forces string method. Figure 2B,C shows the Bayesian measures p(TP|r) for the open-to-semiclosed and semiclosed-to-closed transitions, respectively. Compared with the free-energy profile, these measures look more sensitive to PC dimensions in stringmethod simulations. In Figure 2B, p(TP|r) on the MFEP with the 2D PC lacks sharpness compared with those with the higher dimensional PCs, and a broad peak is wrongly located at around the 22th image. This broadening is explained by the geometric properties of the 2D PC; as shown in Figure 1B, 2D PC cannot clearly divide the semiclosed state from the others. Interestingly, the MFEP with 3D or 4D PCs still has broad distributions, and a sharp peak gradually forms from the MFEP with 10D PC. The peak of p(TP|r) with the 25D PC is sharply located at the 25th image, suggesting that a clear dividing surface with infrequent recrossings is located at this point. For the semiclosed-to-closed transition (Figure 2C), the distribution of the 2D PC has higher values compared with the others, but its shape is too broad to identify its peak position. In contrast, the probability of the 25D PC peaks sharply at the 41th image. Overall, this result indicates that the quality of the MFEP for describing the transition states systematically improves as the dimensionality of PC space increases, 1448

DOI: 10.1021/acs.jpclett.6b00317 J. Phys. Chem. Lett. 2016, 7, 1446−1451

Letter

The Journal of Physical Chemistry Letters

motions of multiple PDB structures while intradomain interactions are kept constant. Here we used the structures of PDB ID 1ake39 and 4ake40 for determining the parameters. Native contact pairs were defined for the open form of AdK (4ake). The effect of ligand binding was incorporated by adding perturbative interactions to the original potential.28 MD simulations were performed with GENESIS,41 a recently developed MD simulation software package in RIKEN AICS. For the PCA, we first performed a brute-force simulation of an AdK DoME model. Using the time-step of 20 fs with SHAKE constraints42 for all bonds, 1.1 × 1010 time steps (corresponding to 220 μs) were integrated. All native contact interactions were calculated without truncation, whereas repulsive nonnative interactions were truncated at a distance of 20 Å. Langevin thermostat was applied for the isothermal condition with a friction coefficient of 0.001 ps−1. Target temperature was 200 K, which is slightly lower than the folding temperature of this model. Outputs were written every 20 ps, resulting in a set of 1.1 × 107 snapshots. During the simulation, 19 924 transitions between open and semiclosed and 1445 transitions between semiclosed and closed were observed. The PCA was performed for the Cartesian coordinates of Cα atoms. Results of the PC modes were used for subsequent string method calculation. The algorithm of mean forces type string method6 was implemented in GENESIS and applied to the open-to-close transition of the AdK DoME model. Because we used the PCs of Cartesian coordinates, tumbling motions were removed by superimposing reference coordinates to the simulation system.26,43 To avoid kinetic trapping, the initial pathway was prepared by piecewise linear lines connecting the average structures of closed, semiclosed, and open states. 64 images were used for a discrete pathway, and the terminal images (first and 64th images) were fixed during the string method simulation. The free-energy profiles were evaluated by projecting the trajectories of the brute-force MD simulation onto the progress coordinate r(x) (eq 1) and calculating the equilibrium densities peq(r). To evaluate samples only near the pathway (MFEP) (i.e., to focus on the samples within a reactive tube6), we introduced a distance from the closest image of the MFEP, which is defined as d(x) = −λ−1 ln-

suggesting that there are risks in characterizing transition-state structures in low-dimensional PC spaces. Then, what kinds of motions exist with higher dimensional PCs and why can the transition states be captured with infrequent recrossings? Figure 2D illustrates structures at the transition state (at the 41th image) between the semiclosed and closed states constructed with the 3D (shown in light colors) and 25D (dark colors) PCs, respectively. A significant difference in the structures is seen around the NMP domain. In the 25D PC, the NMP domain is slightly deformed by intradomain bending, and part of the NMP domain starts to make contact with the LID domain. In contrast, the 3D PC shows no such deformation. Because of this bending within the NMP domain in the 25D PC, the contacts (indicated by lines in the Figure) between the LID and NMP domains form cooperatively and quickly stabilize the closing of the NMP domain and vice versa, resulting in a dividing surface with infrequent recrossings. Interestingly, the same type of NMP domain bending was also observed in the 10 μs trajectory of allatom model28 and the domain detection analysis using two crystal structures,32 suggesting that it is an inherent feature of AdK. We have systematically investigated how the quality of the MFEP calculated by the string method changes by increasing the number of PCs for CV space. For the conformational transitions of AdK, we found that at least 10D PCs are required as PCs to identify the transition states. Comparisons of transition-state structures show that both the global rigid-body domain motion and local intradomain bending are important for capturing a reliable dividing surface. Because local bending motions generally appear in relatively large amplitude PCs, multidimensional PC space should be used with the case of the string method. This is practical because calculation cost does not scale with number of CVs in the string method. Although we anticipate that similar domain deformations would be also observed in other multidomain proteins, transferability of this notion (i.e., importance of domain deformations described by practical number of PCs) to other multidomain proteins will be tested in our future studies. CV choices may include other dimensional reduction techniques, such as time-structure-based independent component analysis (tICA) or relaxation mode analysis.33,34 Local coordinates for each domain by quaternion representation35 could be promising for capturing intradomain bending motions, but in the case of quaternion representation, local domains need to be identified before simulation.32,36 Systematic characterization within these advanced reduction techniques, as performed in this study, may reveal further insights into conformational changes of biomolecules. Finally, the limitation of a coarse-grained model used in this work is noted. In the case of all-atom models, side-chain packings sometimes play important roles on conformational changes. Incorporating side chains in PCs is a challenging issue because the motions of side chains are highly nonlinear and following different metrics from those of backbone. Developments of appropriate transformations37 or kernel type methods38 would be required.

2

−λ∥z(x)−z(m)∥ 30 (∑M ). Here λ is the same as that of eq 1, m=1e the inverse of mean square displacement between successive images. Using this distance from the MFEP, samples that deviated further than a cutoff radius of 10 000 Å2 were ignored. Changing the cutoff radius to 5000 or 15 000 Å2 did not alter the results qualitatively. For the evaluation of the Bayesian measure p(TP|r), the trajectory of the brute-force MD simulation was used. First, trajectory segments in the open, semiclosed, and closed states were assigned according to their RMSD values. Let us define that Ropen is RMSD from the open form and Rclosed is RMSD from the closed form. We assigned the open state when (Ropen − Rclosed) ≤ 3.5 Å, the semiclosed state when 0.5 Å ≤ (Ropen − Rclosed) ≤ 0.5 Å, and the open state when 4.0 Å ≤ (Ropen − Rclosed). Using eq 1, the trajectory was projected onto the progress coordinate for the MFEP, and the probability distributions of r were evaluated. Finally, by applying the Bayesian rule (eq 2), the Bayesian measure p(TP|r) was calculated.

■

COMPUTATIONAL METHODS For MD simulation of AdK, we used a Cα-based coarse-grained model28 recently developed by the authors, called the DoME model. This model enables a stable and efficient MD simulation of multidomain proteins. In the DoME model, the parameters for interdomain interactions are defined by analysis of domain 1449

DOI: 10.1021/acs.jpclett.6b00317 J. Phys. Chem. Lett. 2016, 7, 1446−1451

Letter

The Journal of Physical Chemistry Letters

■

Pathways in Rare Events Simulations. Chem. Phys. Lett. 2006, 426, 168−175. (11) Maragliano, L.; Roux, B.; Vanden-Eijnden, E. Comparison between Mean Forces and Swarms-of-Trajectories String Methods. J. Chem. Theory Comput. 2014, 10, 524−533. (12) Pan, A. C.; Weinreich, T. M.; Shan, Y.; Scarpazza, D. P.; Shaw, D. E. Assessing the Accuracy of Two Enhanced Sampling Methods Using Egfr Kinase Transition Pathways: The Influence of Collective Variable Choice. J. Chem. Theory Comput. 2014, 10, 2860−2865. (13) Best, R. B.; Hummer, G. Reaction Coordinates and Rates from Transition Paths. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 6732−6737. (14) Piana, S.; Lindorff-Larsen, K.; Shaw, D. E. Atomic-Level Description of Ubiquitin Folding. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 5915−5920. (15) Best, R. B.; Hummer, G.; Eaton, W. A. Native Contacts Determine Protein Folding Mechanisms in Atomistic Simulations. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 17874−17879. (16) Abrams, C. F.; Vanden-Eijnden, E. Large-Scale Conformational Sampling of Proteins Using Temperature-Accelerated Molecular Dynamics. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 4961−4966. (17) Vashisth, H.; Brooks, C. L. Conformational Sampling of Maltose-Transporter Components in Cartesian Collective Variables Is Governed by the Low-Frequency Normal Modes. J. Phys. Chem. Lett. 2012, 3, 3379−3384. (18) Jolliffe, I. Principal Component Analysis; Wiley Online Library: 2002. (19) Kitao, A.; Hayward, S.; Go, N. Energy Landscape of a Native Protein: Jumping-among-Minima Model. Proteins: Struct., Funct., Genet. 1998, 33, 496−517. (20) Wolf-Watz, M.; Thai, V.; Henzler-Wildman, K.; Hadjipavlou, G.; Eisenmesser, E. Z.; Kern, D. Linkage between Dynamics and Catalysis in a Thermophilic-Mesophilic Enzyme Pair. Nat. Struct. Mol. Biol. 2004, 11, 945−949. (21) Henzler-Wildman, K. A.; Lei, M.; Thai, V.; Kerns, S. J.; Karplus, M.; Kern, D. A Hierarchy of Timescales in Protein Dynamics Is Linked to Enzyme Catalysis. Nature 2007, 450, 913−916. (22) Hanson, J. A.; Duderstadt, K.; Watkins, L. P.; Bhattacharyya, S.; Brokaw, J.; Chu, J.-W.; Yang, H. Illuminating the Mechanistic Roles of Enzyme Conformational Dynamics. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 18055−18060. (23) Arora, K.; Brooks, C. L. Large-Scale Allosteric Conformational Transitions of Adenylate Kinase Appear to Involve a Population-Shift Mechanism. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 18496−18501. (24) Whitford, P. C.; Miyashita, O.; Levy, Y.; Onuchic, J. N. Conformational Transitions of Adenylate Kinase: Switching by Cracking. J. Mol. Biol. 2007, 366, 1661−1671. (25) Daily, M. D.; Phillips, G. N., Jr.; Cui, Q. Many Local Motions Cooperate to Produce the Adenylate Kinase Conformational Transition. J. Mol. Biol. 2010, 400, 618−631. (26) Matsunaga, Y.; Fujisaki, H.; Terada, T.; Furuta, T.; Moritsugu, K.; Kidera, A. Minimum Free Energy Path of Ligand-Induced Transition in Adenylate Kinase. PLoS Comput. Biol. 2012, 8, e1002555. (27) Lee, J.; Joo, K.; Brooks, B. R.; Lee, J. The Atomistic Mechanism of Conformational Transition of Adenylate Kinase Investigated by Lorentzian Structure-Based Potential. J. Chem. Theory Comput. 2015, 11, 3211−3224. (28) Kobayashi, C.; Matsunaga, Y.; Koike, R.; Ota, M.; Sugita, Y. Domain Motion Enhanced (Dome) Model for Efficient Conformational Sampling of Multidomain Proteins. J. Phys. Chem. B 2015, 119, 14584−14593. (29) Hummer, G. From Transition Paths to Transition States and Rate Coefficients. J. Chem. Phys. 2004, 120, 516−523. (30) Branduardi, D.; Gervasio, F. L.; Parrinello, M. From A to B in Free Energy Space. J. Chem. Phys. 2007, 126, 054103. (31) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009, 180, 1961−1972.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address ⊥

Y.K.: Eisai Co., Ltd., Tsukuba Research Laboratories, Tsukuba-shi, Ibaraki 300-26, Japan.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Computational resources were provided by HOKUSAI GreatWave in RIKEN, RICC (RIKEN Integrated Cluster of Clusters) and K computer in RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Project ID: ra000009). Part of this research has been funded by strategic programs for innovation research, “Computational Life Science and Application in Drug Discovery and Medical Development” (to Y.S.), “Novel Measurement Techniques for Visualizing Live Protein Molecules at Work” (Grant No. 26119006) (to Y.S.), and “Initiative for High-Dimensional Data-Driven Science through Deepening of Sparse Modeling” of MEXT Japan (Grant No. 26120533) (to Y.M.).

■

REFERENCES

(1) Henzler-Wildman, K. A.; Thai, V.; Lei, M.; Ott, M.; Wolf-Watz, M.; Fenn, T.; Pozharski, E.; Wilson, M. A.; Petsko, G. A.; Karplus, M.; Hübner, C. G.; Kern, D. Intrinsic Motions Along an Enzymatic Reaction Trajectory. Nature 2007, 450, 838−844. (2) Shaw, D. E.; Maragakis, P.; Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Eastwood, M. P.; Bank, J. A.; Jumper, J. M.; Salmon, J. K.; Shan, Y.; Wriggers, W. Atomic-Level Characterization of the Structural Dynamics of Proteins. Science 2010, 330, 341−346. (3) Shaw, D. E.; Grossman, J. P.; Bank, J. A.; Batson, B.; Butts, J. A.; Chao, J. C.; Deneroff, M. M.; Dror, R. O.; Even, A.; Fenton, C. H.; Forte, A.; Gagliardo, J.; Gill, G.; Greskamp, B.; Ho, C. R.; Ierardi, D. J.; Iserovich, L.; Kuskin, J. S.; Larson, R. H.; Layman, T.; Lee, L.-S.; Lerer, A. K.; Li, C.; Killebrew, D.; Mackenzie, K. M.; Mok, S. Y.-H.; Moraes, M. A.; Mueller, R.; Nociolo, L. J.; Peticolas, J. L.; Quan, T.; Ramot, D.; Salmon, J. K.; Scarpazza, D. P.; Schafer, U. B.; Siddique, N.; Snyder, C. W.; Spengler, J.; Tang, P. T. P.; Theobald, M.; Toma, H.; Towles, B.; Vitale, B.; Wang, S. C.; Young, C., Anton 2: Raising the Bar for Performance and Programmability in a Special-Purpose Molecular Dynamics Supercomputer. In Proceedings of the 2014 International Conference for High Performance Computing, Networking, Storage and Analysis; IEEE: 2014; pp 41−53. (4) Ohmura, I.; Morimoto, G.; Ohno, Y.; Hasegawa, A.; Taiji, M. MDGRAPE-4: A Special-Purpose Computer System for Molecular Dynamics Simulations. Philos. Trans. R. Soc., A 2014, 372, 20130387. (5) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (6) Maragliano, L.; Fischer, A.; Vanden-Eijnden, E.; Ciccotti, G. String Method in Collective Variables: Minimum Free Energy Paths and Isocommittor Surfaces. J. Chem. Phys. 2006, 125, 024106. (7) Maragliano, L.; Vanden-Eijnden, E. On-the-Fly String Method for Minimum Free Energy Paths Calculation. Chem. Phys. Lett. 2007, 446, 182−190. (8) Pan, A. C.; Sezer, D.; Roux, B. Finding Transition Pathways Using the String Method with Swarms of Trajectories. J. Phys. Chem. B 2008, 112, 3432−3440. (9) Zhang, B. W.; Jasnow, D.; Zuckerman, D. M. The “Weighted Ensemble” Path Sampling Method Is Statistically Exact for a Broad Class of Stochastic Processes and Binning Procedures. J. Chem. Phys. 2010, 132, 054107. (10) Maragliano, L.; Vanden-Eijnden, E. A Temperature Accelerated Method for Sampling Free Energy and Determining Reaction 1450

DOI: 10.1021/acs.jpclett.6b00317 J. Phys. Chem. Lett. 2016, 7, 1446−1451

Letter

The Journal of Physical Chemistry Letters (32) Koike, R.; Ota, M.; Kidera, A. Hierarchical Description and Extensive Classification of Protein Structural Changes by Motion Tree. J. Mol. Biol. 2014, 426, 752−762. (33) Naritomi, Y.; Fuchigami, S. Slow Dynamics in Protein Fluctuations Revealed by Time-Structure Based Independent Component Analysis: The Case of Domain Motions. J. Chem. Phys. 2011, 134, 065101. (34) Mitsutake, A.; Iijima, H.; Takano, H. Relaxation Mode Analysis of a Peptide System: Comparison with Principal Component Analysis. J. Chem. Phys. 2011, 135, 164102. (35) Moradi, M.; Tajkhorshid, E. Computational Recipe for Efficient Description of Large-Scale Conformational Changes in Biomolecular Systems. J. Chem. Theory Comput. 2014, 10, 2866−2880. (36) Hayward, S.; Kitao, A.; Berendsen, H. J. C. Model-Free Methods of Analyzing Domain Motions in Proteins from Simulation: A Comparison of Normal Mode Analysis and Molecular Dynamics Simulation of Lysozyme. Proteins: Struct., Funct., Genet. 1997, 27, 425− 437. (37) Mu, Y.; Nguyen, P. H.; Stock, G. Energy Landscape of a Small Peptide Revealed by Dihedral Angle Principal Component Analysis. Proteins: Struct., Funct., Genet. 2005, 58, 45−52. (38) Schwantes, C. R.; Pande, V. S. Modeling Molecular Kinetics with Tica and the Kernel Trick. J. Chem. Theory Comput. 2015, 11, 600−608. (39) Muller, C. W.; Schulz, G. E. Structure of the Complex between Adenylate Kinase from Escherichia Coli and the Inhibitor Ap5a Refined at 1.9 a Resolution. A Model for a Catalytic Transition State. J. Mol. Biol. 1992, 224, 159−177. (40) Muller, C. W.; Schlauderer, G. J.; Reinstein, J.; Schulz, G. E. Adenylate Kinase Motions During Catalysis: An Energetic Counterweight Balancing Substrate Binding. Structure 1996, 4, 147−156. (41) Jung, J.; Mori, T.; Kobayashi, C.; Matsunaga, Y.; Yoda, T.; Feig, M.; Sugita, Y. GENESIS: A Hybrid-Parallel and Multi-Scale Molecular Dynamics Simulator with Enhanced Sampling Algorithms for Biomolecular and Cellular Simulations. WIREs Comput. Mol. Sci. 2015, 5, 310−323. (42) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. NumericalIntegration of Cartesian Equations of Motion of a System with Constraints - Molecular-Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (43) Branduardi, D.; Faraldo-Gomez, J. D. String Method for Calculation of Minimum Free-Energy Paths in Cartesian Space in Freely-Tumbling Systems. J. Chem. Theory Comput. 2013, 9, 4140− 4154.

1451

DOI: 10.1021/acs.jpclett.6b00317 J. Phys. Chem. Lett. 2016, 7, 1446−1451