Domain Motion Enhanced (DoME) Model for ... - ACS Publications

Nov 4, 2015 - conformational changes between two structures in a protein, provides information on rigid structural units (domains) and the magnitudes ...
0 downloads 8 Views 5MB Size
Download PDF
Article pubs.acs.org/JPCB

Domain Motion Enhanced (DoME) Model for Efficient Conformational Sampling of Multidomain Proteins Chigusa Kobayashi,† Yasuhiro Matsunaga,† Ryotaro Koike,‡ Motonori Ota,‡ and Yuji Sugita*,†,§,∥,⊥ †

RIKEN Advanced Institute for Computational Science, 6-7-1 minatojima-minamimachi, Chuo-ku, Kobe, Hyogo 640-0047, Japan Graduate School of Information Science, Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8601, Japan § RIKEN Theoretical Molecular Science Laboratory and ∥RIKEN iTHES, 2-1 Hirosawa, Wako-shi, Saitama 351-0198, Japan ⊥ RIKEN Quantitative Biology Center, 7-1-26 minatojima-minamimachi, Chuo-ku, Kobe, Hyogo 640-0047, Japan ‡

S Supporting Information *

ABSTRACT: Large conformational changes of multidomain proteins are difficult to simulate using all-atom molecular dynamics (MD) due to the slow time scale. We show that a simple modification of the structure-based coarse-grained (CG) model enables a stable and efficient MD simulation of those proteins. “Motion Tree”, a tree diagram that describes conformational changes between two structures in a protein, provides information on rigid structural units (domains) and the magnitudes of domain motions. In our new CG model, which we call the DoME (domain motion enhanced) model, interdomain interactions are defined as being inversely proportional to the magnitude of the domain motions in the diagram, whereas intradomain interactions are kept constant. We applied the DoME model in combination with the Go model to simulations of adenylate kinase (AdK). The results of the DoME−Go simulation are consistent with an all-atom MD simulation for 10 μs as well as known experimental data. Unlike the conventional Go model, the DoME−Go model yields stable simulation trajectories against temperature changes and conformational transitions are easily sampled despite domain rigidity. Evidently, identification of domains and their interfaces is useful approach for CG modeling of multidomain proteins.

■

INTRODUCTION Proteins exhibit conformational fluctuations around a native structure or transitions between multiple conformers to function under physiological conditions.1,2 X-ray crystallography often provides more than one structure that mimic intermediate states in a reaction cycle.3,4 Large conformational motions of proteins can originate as domain motions, where rigid structural units (domains) change their positions and/or orientations with respect to each other through flexible hinges or loops.5,6 Several computational methods identify multiple domains in a protein from two or more experimental structures7−14 or from simulation trajectories.15−18 A major problem is the definition of domains often depends on how to superimpose one structure on another, for instance, which regions of two structures are superimposed with each other. Recently, a hierarchical structural comparison method, which is free from the superimposition problem, has been developed by Koike et al.19 A tree diagram called “Motion Tree”19 describes conformational changes of a protein in a hierarchical manner. Molecular dynamics (MD) simulations of proteins have been widely used to investigate conformational dynamics, stability and structure−function relationships. Either an all-atom model20−25 or a coarse-grained (CG) model is utilized.26−35 MD simulations based on an all-atom model provide the most accurate descriptions of protein structures and dynamics. © XXXX American Chemical Society

Simulation runs are usually not longer than microseconds, which may not be sufficient to observe domain motions in a protein. In contrast, MD simulations based on a CG model extend the time scale to milliseconds, which covers longer time scale motions, such as protein folding,36−39 aggregation,40−42 and large-scale conformational changes in proteins or protein complexes.43−46 CG MD simulations (and their thermodynamics and kinetics) are strongly dependent on CG force-field parameters and their selection is critical. A simple off-lattice Go-like potential, which is based on the funnel theory of protein folding and dynamics,47−49 has often been used in CG MD simulations.28,50 Here, interactions between residues in close proximity in a native structure, called “native contacts”, are the important nonbonded interactions. More recently, the native contact parameters have been modified to improve the accuracy of CG MD simulations.29−31,33,34 Karanicolas and Brooks adapted knowledgebased parameters (the Miyazawa−Jernigan contact energy51) to the native contacts and succeeded in distinguishing alternate folding pathways of different proteins with the same topology.29 Hereafter, we refer to it as the KB−Go model. Other recent Received: August 7, 2015 Revised: September 28, 2015

A

DOI: 10.1021/acs.jpcb.5b07668 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B approaches are force-matching,30,31 fluctuation-matching,34 or REACH33 methods, where all-atom MD simulations are performed in advance and the forces, fluctuations, or covariances of a CG model are fitted to those of all-atom simulations. These algorithms allow simulation of not only protein folding/unfolding but also many other essential biological phenomena.45,46,52−54 Here, we propose a new CG model for simulating large-scale motions in multidomain proteins. First, domains and magnitudes of domain motions in a protein are identified in the Motion Tree.19 Then, crucially, the interdomain nativecontact parameters are defined as being inversely proportional to the magnitudes of domain motions in the diagram. This CG modeling, which we call DoME (domain motion enhanced) model, enhances conformational sampling of multidomain proteins. We applied our DoME model combined with the Go model (DoME−Go model) to a small soluble globular protein, adenylate kinase (AdK). AdK has 214 residues and consists of three domains (LID, NMP, and CORE). The protein catalyzes the following reaction,

distance in the native structure and the minimum depths of the interaction, εij, are determined empirically. In the KB−Go model, εij is given as the Miyazawa−Jernigan contact energy51 for each residue type, which was derived by statistical analysis of PDB structures. The reference parameters of the dihedral angle term are also given as knowledge-based ones in the KB−Go model. Dual-Basin Go Model Based on the Perturbation Approach. The Go model described in eq 1 is categorized as a single-basin model, where the native contacts in the potential energy function are derived from a single native structure. It is especially useful for protein folding/unfolding simulations, whereas it may have difficulty in simulating transitions between two states. In contrast, two or more structures are stabilized in the potential energy function in multibasin Go models.32,35,61,69−74 Whitford et al. developed a perturbation approach61 where a single-basin Go model is perturbed by additional native contacts that exist uniquely in the other structure. We here utilize this approach to study transitions between the open and closed forms of AdK and other multidomain proteins. Starting from a single-basin Go model for the open form, the perturbation toward the closed form is then added. The potential energy function for the perturbation is equivalent to that of the native contacts in eq 1. We performed several preliminary simulations by changing the depth of the perturbation interaction, εLigand, and then determined the optimal parameters for production dynamics. Domain Motion Enhanced (DoME) Model for Multidomain Proteins. We developed the DoME model for multidomain proteins that have more than one experimental structure. This model utilizes the Motion Tree that describes conformational changes in a protein, and we provide a brief account of it here. Two different structures (a and b) of a protein consisting of N residues are considered. We compute the difference of two distance matrices eij = |daij − dbij|, where dXij (i, j = 1, ..., N and X = a, b) is the Cα-atom distance matrix for each structure. A tree diagram named Motion Tree is computed by hierarchical clustering of eij. Internal nodes in the tree represent rigid-body domains in various scale-size ranges. By construction, each ancestral node contains two descendant nodes, which can be regarded as smaller domains responsible for internal motions within a larger domain. The average distance difference between these two domains represents the magnitude of rigid-body motion, denoted by M hereafter (see precise definition in ref 19). The M value is indicated by the tree height at the ancestral node of two domains. Following the original study, we consider a magnitude of 5.0 Å19 as the criteria to decide whether displacements are considerable as rigid-body domain motions or not. In other words, we regard two residues as being within the same rigid body if the height at their common ancestor (the magnitude) is less than 5.0 Å and discard any internal motions in the rigid body (see Figure 1A in detail). Hereafter, we use the rigid domain, or just the domain as the indivisible rigid body in this sense above. The DoME model combined with the Go model (DoME− Go) has the same potential function of the KB−Go model except for the nonbonded native contact parameters, εij, between atoms i and j. We set the constant value of εij for the intradomain atom pairs and modify εij only for the interdomain atom pairs, as shown in eq 2.

ATP + AMP ⇄ 2ADP

AdK adopts a closed form with bound ligands and an open structure in the apo state.55 This conformational transition is the rate-limiting step in the enzymatic reaction.3,56 The mechanism of open-to-closed conformational transition has been extensively investigated experimentally3,55−60 and in theoretical/computational studies.61−66 We compared the fluctuations of DoME−Go with those of KB−Go and an allatom explicit solvent MD trajectory for 10 μs using a GPGPU computer.67,68 We also simulated the conformational transition between the open and closed forms in AdK and other multidomain proteins.61 We indeed observe many conformational transitions in the DoME−Go simulations, suggesting the usefulness of our model in structural modeling of multidomain proteins.

■

COMPUTATIONAL METHODS Potential Energy Function of the Off-Lattice Go Model. The potential energy function of the off-lattice Go model consists of bonded terms as well as nonbonded terms: E=

∑ Kb(b − b0)2 + ∑ Kθ(θ − θ)2 bond

+

∑

angle

Kϕ(1 + cos(nϕ − ϕ0))

dihedral

⎡ ⎛ ⎞12 ⎛ σij ⎞10 ⎛ σij ⎞6 ⎤ σij ⎢ ⎜ ⎟ ⎜ ⎟ ∑ εij⎢13⎜ ⎟ − 18⎜ ⎟ + 4⎜⎜ ⎟⎟ ⎥⎥ + ⎝ rij ⎠ ⎝ rij ⎠ ⎦ native contacts ⎣ ⎝ rij ⎠ ⎡⎛ ⎞12 ⎤ R ij ∑ + ε⎢⎜⎜ ⎟⎟ ⎥ ⎢⎝ rij ⎠ ⎥ non‐native contacts ⎣ ⎦

(1)

The reference parameters b0, θ0, and ϕ0 for the bond, angle, and dihedral terms are taken from the native structure of a target protein. The nonbonded interaction consists of native and nonnative contact interactions. The former is computed for atom pairs whose distances in the native structure are within predefined criteria, whereas the latter is computed for the rest of the pairs as a repulsive interaction. σij is taken from the B

DOI: 10.1021/acs.jpcb.5b07668 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B ⎧ c (i ∈ I , j ∈ J , I ≠ J ) ⎪ ⎪ MIJ εij = ⎨ c ⎪ c = (i , j ∈ I) ⎪M M ⎩ II criterion

(2)

εij for the interdomain atom pairs is defined as inversely proportional to the magnitude, MIJ, where atoms i and j belong to domains I and J, respectively. Mcriterion is set to 5.0 Å, the criteria used in the original study.19 A constant factor c is determined to keep a reasonable balance between the nonbonded native contacts and the rest of the interactions. This model is based on the assumption that the magnitude of interdomain motions is mainly determined by interactions at the domain interface. This assumption is consistent with the hierarchical descriptions of domain motions by Motion Tree. The internal motion of a large part is interpreted as the motions of two smaller domains. This suggests that the interface between the smaller domains is more flexible compared to other regions. The idea of modeling interactions at the domain interfaces is supported by a number of normal-mode analysis studies.8,75,76 For example, Lu and Ma showed that the largeamplitude or low-frequency motions are robust upon randomization of the Hessian matrix elements as long as its block diagonal structure (where each block corresponds to a rigidbody domain) is maintained.76 Simulation Details. We examined the efficiency and stability of MD simulations based on the DoME−Go model, monitoring the fluctuations of AdK in the open form and the transition between the open and closed forms. The former was simulated using a single-basin model with 476 native contact pairs, whereas the latter was simulated with a dual-basin model with 17 additional perturbative contacts toward closed form. We also carried out MD simulations based on the KB−Go model as well as an all-atom explicit solvent MD simulation using AMBER force fields.23,77,78 Two crystal structures of AdK (PDBID: 4AKE55 and 1AKE58) were used for the open and closed forms, respectively, to set up the DoME−Go and KB−Go models. Only Cα atoms were used. All CG MD simulations were performed using GENESIS,79 which has been developed in RIKEN AICS. A Langevin thermostat that realizes isothermal ensembles was used for 10 μs with a friction coefficient of 0.001 ps−1. All bonds were fixed with the SHAKE algorithm,80 and the time step was 0.02 ps. All native contact interactions were calculated without truncation, whereas repulsive non-native interactions were truncated at a distance of 20.0 Å. We evaluated heat capacities using WHAM81 on the trajectories of DoME−Go and KB−Go simulations at different temperatures. AdK showed multiple peaks in heat capacity, which correspond to different folding temperatures of three domains. The existence of multiple peaks has been reported in a previous CG simulation.82 We defined the lowest folding temperature obtained in KB−Go simulation as Tf. We mainly set 0.9Tf in the production runs of DoME−Go and KB−Go simulations. In addition, five more temperatures (1.0Tf, 1.1Tf, 1.2Tf, 1.3Tf, and 1.4Tf) were used to examine the conformational stability of AdK. We also performed an all-atom MD simulation starting from the open structure for 10 μs, using AMBER software package (version 12) on a workstation with GPGPU (GTX TITAN, SPFP precision model68,83). AMBER ff12SB23,77,78 and the TIP3P model were used for the protein and solvent water,

Figure 1. Definition of rigid domains in the Motion Tree (MT) and parameters in the DoME−Go model for adenylate kinase (AdK). (A) MT and rigid domains in the reaction step between the open and closed forms of AdK. The top left cartoon represents the division of rigid domains as clusters at the threshold of 5.0 Å magnitudes (red dotted line in the MT). The bottom protein structures represent the domain separations in AdK. The effective node number and the magnitude in MT are written above the protein structures. Two different colors other than gray mean that AdK is divided into two parts on each effective node. The residues in gray indicate that they are not included in the separation on the node. The top and bottom structures are taken from the crystal structures in closed form with ligand and open form, respectively. (B) Magnitudes of interdomains, MIJ in eq 2. The native contact parameters (εij in eq 1) in the potential functions of (C) the DoME−Go model and (D) the KB−Go model, both mapped onto the contact map of AdK in the open form.

C

DOI: 10.1021/acs.jpcb.5b07668 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

The native contact parameters εij in the DoME−Go model are compared with those of the KB−Go model in Figure 1C,D. Dotted lines in the figures indicate domains defined in the Motion Tree. The Miyazawa−Jernigan parameters used in the KB−Go model have strong interactions between hydrophobic residue pairs, whereas the hydrophilic interactions are weaker. By construction, all the intradomain contact parameters in DoME−Go are constant, and the parameters in the hinge regions connecting to the rigid domains are weak. We defined the constant factor in eq 2 such that the intrahelix interaction in the DoME−Go model is the same as that in the KB−Go model. All-Atom MD Simulation of AdK Starting from the Open Form. In the all-atom MD simulation, the whole structure of AdK remained stable during 10 μs (Figure 2A). We first compare the root-mean-square fluctuation (RMSF) of the Cα atoms for 10 μs trajectory (RMSF10 μs) with that for the initial 10 ns (RMSF10 ns) in Figure 2B. RMSF10 μs of LID1 are about 1.5 times larger than RMSF10 ns in the same regions. The fluctuation of LID1 is smaller than that of NMP for the initial 10 ns, whereas it becomes larger in RMSF10 μs. We examined the time dependency of large-amplitude motions in AdK by performing principal component analysis (PCA) on the trajectories obtained for 10 μs, the initial 1 μs, 100 ns, and 10 ns. In Figure 2C−E, the 20 largest-amplitude PC vectors for 10 μs and those for the initial 1 μs, 100 ns, and 10 ns are compared by taking their inner products. Figure 2C shows a good correspondence between the PC vectors for 10 μs and the initial 1 μs, and Figure 2D,E shows weaker correlations, suggesting that at least 1 μs all-atom MD simulation is necessary to evaluate conformational fluctuations of AdK, properly. Conformational Fluctuations of AdK in the Open Form by DoME−Go Simulations. We evaluated conformational fluctuations and the stability of AdK in simulations based on DoME−Go models and compared them to those in the allatom and KB−Go simulations. To investigate conformational

respectively. The simulation system consists of an AdK, 3 Na+ atoms, and 12 336 water molecules to fill the simulation box of size 68.9 Å × 86.8 Å × 83.8 Å. The particle mesh Ewald method84,85 was applied to long-range electrostatic interactions; the short-range Coulomb and van der Waals interactions are truncated at 10 Å. In a production run, temperature and pressure are controlled at 300 K and 1 atm by the Langevin thermostat and the Berendsen barostat,86 respectively.

■

RESULTS Native Contact Interaction in the DoME−Go Model. In Figure 1A, a Motion Tree diagram of the conformational changes of AdK between open and closed forms and three domain motions whose magnitudes are greater than 5.0 Å is shown. The first and second nodes in the diagram describe the motions involving head loops of the LID and of the NMP, respectively. The magnitudes of the two nodes are 23.3 and 19.7 Å. The third node has a smaller magnitude (8.4 Å) and indicates the motions of two helices in the LID and P-loop (residues 7−13) relative to the CORE. In this case, 4 domains (LID loops, NMP, LID helices, and rest) were identified (Figure 1A and Table 1). The MIJ values employed in eq 2 for Table 1. Residues in Separated Domains of Adenlylate Kinase Defined by Motion Tree domain

residues

CORE LID1 LID2 NMP

1−6, 14−29, 60−114, 175−214 121−158 7−13, 115−120, 159−174 30−59

the intradomain interaction are shown in Figure 1B. The P-loop is a functionally important motif, which specifically binds ATP. We hereinafter refer to the head loops of the LID as LID1, whereas the two helices and P-loop are referred to as LID2.

Figure 2. Conformational fluctuations from the all-atom MD simulation trajectory. (A) RMSD from the crystal structure in the open form as a function of time. (B) RMSF from the trajectory for 10 μs (solid line) and the initial 10 ns (dotted line), respectively. (C)−(E) Inner products of the 20 lowest frequency principal component (PC) vectors in 10 μs trajectory and those in (C) the initial 1 μs, (D) the initial 100 ns, and (E) the initial 10 ns, respectively. D

DOI: 10.1021/acs.jpcb.5b07668 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 3. Conformational fluctuations and stabilities of the AdK structure in the all-atom, DoME−Go, and KB−Go simulations. RMSFs (A)−(C), fractions of the native contacts (Q-value) of a whole AdK structure (D)−(F), and LID1 (G)−(I) are shown. The top, middle, and bottom figures show the results of the all-atom MD (A, D, and G), DoME−Go (B, E, and H), and KB−Go (C, F, and I) simulations. Black, blue, green, yellow, orange, and red in (B, C, E, F, H, and I) represent results obtained at different temperatures, 0.9Tf, 1.0Tf, 1.1Tf, 1.2Tf, 1.3Tf, and 1.4Tf, respectively.

Figure 4. Variances of distances between Cα atoms in the (A) all-atom, (B) DoME−Go, and (C) KB−Go simulations.

diagonal regions and large variances in the distance between the LID1 and NMP domains. In contrast, the variances of the nondiagonal regions in the KB−Go simulations are small. The average of the diagonal regions (residues 121−158) in the KB− Go simulation is 1.09 Å2, much greater than those in the allatom MD (0.35 Å2) and DoME−Go simulations (0.26 Å2) . In summary, DoME−Go and all-atom MD simulations provided consistent results on the conformational fluctuations in AdK. The partial unfolding of the LID1 domain was observed only in the KB−Go model, presumably due to the weaker hydrophilic interactions in the model (Figure 1D). Conformational Transition of AdK between the Open and Closed Forms. We investigated the conformational transition of AdK between the open and closed forms using the dual-basin Go model in combination with DoME and KB models (see Figure S1 for εij). We used the two angles, θNMP and θLID, as reaction coordinates to draw the free-energy landscape of the conformational transition in Figures 5 and S2. θNMP is the angle defined by the three centers of mass for the residues 115−125, 90−100, and 35−55. θLID is the angle by the centers of mass for the residues 179−185, 115−125, and 125−

stability of the models against temperature changes, we performed DoME−Go and KB−Go simulations at six different temperatures (0.9Tf, 1.0Tf, 1.1Tf, 1.2Tf, 1.3Tf, and 1.4Tf). In Figure 3, RMSFs and fractions of native contacts, so-called Qvalues,28 are shown for the all-atom MD, DoME−Go, and KB− Go simulations. If large RMSFs are observed in a protein in spite of Q-values being close to 1, the motion of the whole protein is mainly attributed to rigid-body domain motions. Large fluctuations in the LID and NMP are observed in DoME−Go simulations at all six temperatures whereas Qvalues in all the simulations are close to 1. This suggests that AdK undergoes domain motions regardless of temperature changes. In contrast, KB−Go simulations at 1.0Tf or higher show large fluctuations in the LID1 domain. At temperatures higher than 1.1Tf, the Q-value of the whole structure is about 0.8 and that of the LID1 domain below 0.5. In Figure 4, the variances of the distances between Cα atoms are shown for the all-atom MD simulation at 300 K, the DoME−Go simulation at 0.9Tf, and the KB−Go simulation at 0.9Tf. The results of DoME−Go simulations are similar to those of all-atom MD simulations: small variances in the E

DOI: 10.1021/acs.jpcb.5b07668 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 6. Potential mean force (PMF) as a function of difference of RMSDs from the open and closed structures with different models (ΔRMSD = RMSDopen − RMSDclosed). (A) PMFs obtained in all-atom (dashed line), DoME−Go (thick line), and KB−Go (thin line) simulations. (B) PMFs obtained in dual-basin Go simulations. Thick and thin lines indicate the results of dual-basin DoME−Go and dualbasin KB−Go simulations, respectively.

The semi-open form, which is only observed in the DoME−Go simulation, is consistent with previous simulation results.32,34,61−65 We further examine the perturbation parameter dependence in Figure 7. Even with a rather small value for the perturbation, DoME−Go simulations allowed 10 transitions between the open and closed forms, whereas during KB−Go simulations the protein remained in the open form. We therefore increased the perturbation interaction (εLigand) systematically to lower the free energy of the closed form. Due to the high energy barrier between the two forms in the KB−Go model, however, conformational transitions between the open and semi-open forms were observed only with εLigand = 2.1 and 2.4, whereas only the semi-open form was stabilized with εLigand = 2.7. With εLigand = 3.2, the closed form becomes dominant, underlying the difficulty in achieving the conformational transition in KB−Go simulations.

Figure 5. Free-energy surfaces of conformational transitions in AdK in (A) the all-atom MD, (B) DoME−Go, and (C) dual-basin DoME−Go simulations. Vertical and horizontal axes indicate the orientation of LID and NMP with respect to CORE, respectively. The definitions are written in text. Red and magenta dots represent the crystal structures in the open and closed forms, respectively. (D) Structures of AdK in the three centroids on the free energy surface by dual-basin DoME− Go simulation. (White triangles in (C)) Open, semi-open, and closed structures are displayed.

■

DISCUSSION AND CONCLUSIONS Mechanism for Reducing the Energy Barrier of Conformational Change. DoME−Go simulations show that the LID undergoes large amplitude motions. This flexibility is consistent with the results of NMR3 and single molecule FRET.59,60 The mechanism for reducing the energy barrier between the open and closed forms in AdK has been discussed in several CG MD studies: Whitford et al. simulated the transition on the basis of a dual-basin Go model with perturbation. Local unfolding (cracking) at the interfaces between domains (i.e., CORE-LID and CORE-NMP) and the accompanying increase of entropy was proposed for a mechanism to decrease the energy barrier.61 In contrast, Daily et al.63 also simulated AdK using the KB−Go model with the macro-mixing approach70 and concluded that cracking may be an artifact in the dihedral angle potential used in the study of Whitford et al. In our study, we show that weakening the interactions at the domain interfaces helps to decrease the energy barrier. This points to the importance of tuning dihedral and nonbonded interactions together for modeling conformational transitions in proteins. Importance of Heterogeneous Native Contact Parameters Derived from Motion Tree. To investigate the importance of heterogeneous native contact parameters given by a hierarchical representation in the Motion Tree, we performed MD simulations by using constant native contact parameters for all the domain interfaces. We constructed four

153. The definition is similar to that used in previous studies.87,88 In the all-atom MD simulation at 300 K, we observed a wide distribution around the crystal structure in the open form and no transition from the open to closed forms (Figure 5A), suggesting the simulation length was not sufficiently long for the transition. No transition from the open to closed forms occurred in KB−Go simulation even with the perturbation in the dual-basin Go model as shown in Figure S2B. The perturbation toward the closed form does, in fact, make for a wider distribution, but it is not sufficient to realize the structural transition within accessible simulation time. In contrast, the DoME−Go simulation provided three stable distributions, as shown in Figure 5C. The centers of the first two distributions are close to the X-ray crystal structures of open and closed forms. The third takes a conformation in which only the LID1 domain is closed and the NMP domain is open. We refer to this structure as the semi-open form. The one-dimensional potential of mean force (PMF) for the transitions (parts A and B of Figure 6 are singlebasin and dual-basin models, respectively) suggests that the free-energy barrier between the open and the semi-open forms is lower than that between the semi-open and closed forms. F

DOI: 10.1021/acs.jpcb.5b07668 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 7. Time series of difference of RMSDs from the open and closed forms with different strengths of perturbation. (A) DoME−Go model with εLigand = 1.7 kcal/mol. (B)−(F) different parameters of perturbation in the KB−Go model. A label in graph indicates εLigand.

models of AdK with MIJ = 23.3, 19.7, 8.4, and 5.0, which were taken from the magnitudes of nodes in the Motion Tree (Figure 1B). In Figure S3, larger MIJ facilitated the amplitude of domain motions in AdK, whereas smaller MIJ kept AdK more rigid. The free-energy surface in Figure S4 shows an additional intermediate state at ΔRMSD ∼ 3 Å in the simulations with MIJ = 23.3 and destabilization of the closed conformations in comparison with the DoME−Go simulations. This change results from the overestimation of NMP flexibility in this model, which was prohibited in the original DoME−Go simulations using the stronger native contact parameters of the interface between NMP and CORE domains. For MIJ = 19.7, the semi-open form was slightly overstabilized compared to the open form. In the simulations with MIJ = 8.4 and 5.0, the open form was overstabilized due to too strong native contact parameters of all the domain interfaces, as shown in Figure S4. In contrast to these models with constant native contact parameters for the domain interfaces, the simulations based on the DoME with heterogeneous native contact parameters provided results consistent with the previous simulations as well as experimental studies, indicating the importance of a subtle balance in native contact parameters for describing the energetics of conformational change. A reasonable set of the heterogeneous parameters obtained by a hierarchical representation in the Motion Tree is necessary for the modeling of the domain motions. General Applicability of the DoME−Go Model to Other Multiple-Domain Proteins. To examine the general applicability of the DoME model for conformational fluctuations, we applied DoME−Go simulations to two other proteins (T4 lysozyme and lysine-, arginine-, ornithine-binding protein (LAO)). The parameters of the DoME−Go model (with a single-basin model) in T4 lysozyme (Figure S5) and LAO (Figure S6) were generated using the procedure written in the Computational Method section. First, except for the lowest temperature, the KB−Go model shows large RMSFs and small Q-values for both proteins, suggesting that unfolding events take place easily. In contrast, the temperature effect on RMSF and Q-values in DoME−Go simulations is much more minor. In particular, the fluctuations of T4 lysozyme were very small. Motion Trees of the two proteins indicate only one smaller domain motion with respect to the larger domain, if we use the

same criteria for defining large domain motions. This allows simpler conformational motions in CG MD simulations based on the DoME−Go model. Advantages and Limitations of the DoME Model. The DoME model has at least two advantages over other structurebased CG modeling of multidomain proteins. The first one is the simplicity of setting up the model parameters. If two or more experimental structures are known for the target protein, the calculation of the Motion Tree and the setting up the DoME model with a single-basin or a dual-basin Go model is automatic without the input of unreliable parameters. This contrasts with other methods, in particular, the forcematching30,31 and fluctuation-matching34 methods use information from all-atom MD simulations. As we have seen, it requires microsecond or longer to converge the RMSFs and the lowest PC vectors of AdK. Another advantage is the robustness of CG MD simulations against temperature changes. We were able to simulate large-amplitude domain motions of AdK, T4 lysozyme, and LAO at different temperatures in DoME−Go simulations, which was not possible in KB−Go simulations. In CG MD simulations, selection of a proper temperature is one of the key setups.28 The robustness of DoME−Go against temperature changes suggests that it will be straightforward to apply the model to temperature replica-exchange MD (TREMD)89,90 for further enhancement of conformational sampling of multidomain proteins. In the case of large proteins consisting of multiple domains, the combination of T-REMD with the DoME model should help sampling wider conformational space. The DoME model has, like other structure-based CG models, limitations. The method is based on the Motion Tree that requires more than one structure, and some of proteins are not available. Detailed analyses of intradomain motions or folding/unfolding events are beyond the scope of the model. There is no general parameter set or model that is applicable to all biological phenomena, like protein folding/ unfolding, protein binding and aggregation, large-amplitude domain motions, etc., and each problem needs to be assessed individually. For instance, the KB model without modification may not be appropriate for simulating multiple-domain proteins, but it is indeed suitable for simulating protein folding, as shown in the original paper.29 In a future application, we plan G

DOI: 10.1021/acs.jpcb.5b07668 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(8) 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. (9) Wriggers, W.; Schulten, K. Protein Domain Movements: Detection of Rigid Domains and Visualization of Hinges in Comparisons of Atomic Coordinates. Proteins: Struct., Funct., Genet. 1997, 29, 1−14. (10) Hayward, S.; Berendsen, H. J. C. Systematic Analysis of Domain Motions in Proteins from Conformational Change: New Results on Citrate Synthase and T4 Lysozyme. Proteins: Struct., Funct., Genet. 1998, 30, 144−154. (11) Verbitsky, G.; Nussinov, R.; Wolfson, H. Flexible Structural Comparison Allowing Hinge-Bending, Swiveling Motions. Proteins: Struct., Funct., Genet. 1999, 34, 232−254. (12) Shatsky, M.; Nussinov, R.; Wolfson, H. J. Flexible Protein Alignment and Hinge Detection. Proteins: Struct., Funct., Genet. 2002, 48, 242−256. (13) Shatsky, M.; Nussinov, R.; Wolfson, H. J. Multiprot - a Multiple Protein Structural Alignment Algorithm. Lect. Notes. Comput. Sci. 2002, 2452, 235−250. (14) Salem, S.; Zaki, M. J.; Bystroff, C. Flexsnap: Flexible NonSequential Protein Structure Alignment. Algorithms Mol. Biol. 2010, 5, 12. (15) Lange, O. F.; Grubmuller, H. Generalized Correlation for Biomolecular Dynamics. Proteins: Struct., Funct., Genet. 2006, 62, 1053−1061. (16) Hub, J. S.; de Groot, B. L. Detection of Functional Modes in Protein Dynamics. PLoS Comput. Biol. 2009, 5, e1000480. (17) Menor, S. A.; de Graff, A. M. R.; Thorpe, M. F. Hierarchical Plasticity from Pair Distance Fluctuations. Phys. Biol. 2009, 6, 036017. (18) Yesylevskyy, S. O. Identifying the Hierarchy of Dynamic Domains in Proteins Using the Data of Molecular Dynamics Simulations. Protein Pept. Lett. 2010, 17, 507−516. (19) 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. (20) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.; Weiner, P. A New Force-Field for Molecular Mechanical Simulation of Nucleic-Acids and Proteins. J. Am. Chem. Soc. 1984, 106, 765−784. (21) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. An All Atom Force-Field for Simulations of Proteins and Nucleic-Acids. J. Comput. Chem. 1986, 7, 230−252. (22) Jorgensen, W. L.; Tiradorives, J. The Opls Potential Functions for Proteins - Energy Minimizations for Crystals of Cyclic-Peptides and Crambin. J. Am. Chem. Soc. 1988, 110, 1657−1666. (23) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A, 2nd Generation Force-Field for the Simulation of Proteins, Nucleic-Acids, and Organic-Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (24) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (25) Scott, W. R. P.; Hunenberger, P. H.; Tironi, I. G.; Mark, A. E.; Billeter, S. R.; Fennen, J.; Torda, A. E.; Huber, T.; Kruger, P.; van Gunsteren, W. F. The Gromos Biomolecular Simulation Program Package. J. Phys. Chem. A 1999, 103, 3596−3607. (26) Tirion, M. M. Large Amplitude Elastic Motions in Proteins from a Single-Parameter, Atomic Analysis. Phys. Rev. Lett. 1996, 77, 1905− 1908. (27) Bahar, I.; Atilgan, A. R.; Erman, B. Direct Evaluation of Thermal Fluctuations in Proteins Using a Single-Parameter Harmonic Potential. Folding Des. 1997, 2, 173−181. (28) Clementi, C.; Nymeyer, H.; Onuchic, J. N. Topological and Energetic Factors: What Determines the Structural Details of the

to investigate enzymatic reactions by examining large-amplitude and long-time scale conformational changes in a multiresolution approach. Multiresolution simulations, combining the DoME model and all-atom modeling, would smooth over weaknesses in each approach and provide insights into protein conformational changes accompanying chemical reactions.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b07668. Six figures showing the parameters of the dual-basin DoME model, the potential of mean force of the KB−Go model, the conformational fluctuations with different interdomain parameters, the potential of mean force with different interdomain parameters, the domain description and conformational fluctuations of the T4 lysozyme, and the domain description and conformational fluctuations of LAO (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Yuji Sugita. E-mail: [email protected]. 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” (Grand No. 26119006) (to Y.S.), “Initiative for HighDimensional Data-Driven Science through Deepening of Sparse Modeling” of MEXT Japan (Grant No. 26120533) (to Y.M.), and “Platform for Drug Discovery, Informatics, and Structural Life Science” from Japan Agency for Medical Research and development (AMED) (to M.O.).

■

REFERENCES

(1) Boehr, D. D.; Nussinov, R.; Wright, P. E. The Role of Dynamic Conformational Ensembles in Biomolecular Recognition. Nat. Chem. Biol. 2009, 5, 789−796. (2) Ma, B. Y.; Nussinov, R. Enzyme Dynamics Point to Stepwise Conformational Selection in Catalysis. Curr. Opin. Chem. Biol. 2010, 14, 652−659. (3) 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. (4) Henzler-Wildman, K.; Kern, D. Dynamic Personalities of Proteins. Nature 2007, 450, 964−972. (5) Wodak, S. J.; Janin, J. Location of Structural Domains in Proteins. Biochemistry 1981, 20, 6544−6552. (6) Gerstein, M.; Lesk, A. M.; Chothia, C. Structural Mechanisms for Domain Movements in Proteins. Biochemistry 1994, 33, 6739−6749. (7) Boutonnet, N. S.; Rooman, M. J.; Wodak, S. J. AutomaticAnalysis of Protein Conformational-Changes by Multiple Linkage Clustering. J. Mol. Biol. 1995, 253, 633−647. H

DOI: 10.1021/acs.jpcb.5b07668 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Transition State Ensemble and ″En-Route″ Intermediates for Protein Folding? An Investigation for Small Globular Proteins. J. Mol. Biol. 2000, 298, 937−953. (29) Karanicolas, J.; Brooks, C. L. The Origins of Asymmetry in the Folding Transition States of Protein L and Protein G. Protein Sci. 2002, 11, 2351−2361. (30) Izvekov, S.; Parrinello, M.; Burnham, C. J.; Voth, G. A. Effective Force Fields for Condensed Phase Systems from Ab Initio Molecular Dynamics Simulation: A New Method for Force-Matching. J. Chem. Phys. 2004, 120, 10896−10913. (31) Izvekov, S.; Voth, G. A. A Multiscale Coarse-Graining Method for Biomolecular Systems. J. Phys. Chem. B 2005, 109, 2469−2473. (32) Maragakis, P.; Karplus, M. Large Amplitude Conformational Change in Proteins Explored with a Plastic Network Model: Adenylate Kinase. J. Mol. Biol. 2005, 352, 807−822. (33) Moritsugu, K.; Smith, J. C. Coarse-Grained Biomolecular Simulation with Reach: Realistic Extension Algorithm Via Covariance Hessian. Biophys. J. 2007, 93, 3460−3469. (34) Li, W. F.; Wolynes, P. G.; Takada, S. Frustration, Specific Sequence Dependence, and Nonlinearity in Large-Amplitude Fluctuations of Allosteric Proteins. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 3504−3509. (35) Terada, T. P.; Kimura, T.; Sasai, M. Entropic Mechanism of Allosteric Communication in Conformational Transitions of Dihydrofolate Reductase. J. Phys. Chem. B 2013, 117, 12864−12877. (36) Tozzini, V. Coarse-Grained Models for Proteins. Curr. Opin. Struct. Biol. 2005, 15, 144−150. (37) Clementi, C. Coarse-Grained Models of Protein Folding: Toy Models or Predictive Tools? Curr. Opin. Struct. Biol. 2008, 18, 10−15. (38) Hyeon, C.; Thirumalai, D. Capturing the Essence of Folding and Functions of Biomolecules Using Coarse-Grained Models. Nat. Commun. 2011, 2, 487. (39) Kamerlin, S. C. L.; Vicatos, S.; Dryga, A.; Warshel, A. CoarseGrained (Multiscale) Simulations in Studies of Biophysical and Chemical Systems. Annu. Rev. Phys. Chem. 2011, 62, 41−64. (40) Nguyen, P. H.; Li, M. S.; Stock, G.; Straub, J. E.; Thirumalai, D. Monomer Adds to Preformed Structured Oligomers of a Beta-Peptides by a Two-Stage Dock-Lock Mechanism. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 111−116. (41) Pellarin, R.; Guarnera, E.; Caflisch, A. Pathways and Intermediates of Amyloid Fibril Formation. J. Mol. Biol. 2007, 374, 917−924. (42) Bellesia, G.; Shea, J. E. Effect of Beta-Sheet Propensity on Peptide Aggregation. J. Chem. Phys. 2009, 130, 145103. (43) Sanbonmatsu, K. Y.; Joseph, S.; Tung, C. S. Simulating Movement of Trna into the Ribosome During Decoding. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 15854−15859. (44) Koga, N.; Takada, S. Folding-Based Molecular Simulations Reveal Mechanisms of the Rotary Motor F-1-Atpase. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 5367−5372. (45) Yao, X. Q.; Kenzaki, H.; Murakami, S.; Takada, S. Drug Export and Allosteric Coupling in a Multidrug Transporter Revealed by Molecular Simulations. Nat. Commun. 2010, 1, 117. (46) Yao, X. Q.; Kimura, N.; Murakami, S.; Takada, S. Drug Uptake Pathways of Multidrug Transporter Acrb Studied by Molecular Simulations and Site-Directed Mutagenesis Experiments. J. Am. Chem. Soc. 2013, 135, 7474−7485. (47) Bryngelson, J. D.; Onuchic, J. N.; Socci, N. D.; Wolynes, P. G. Funnels, Pathways, and the Energy Landscape of Protein-Folding - a Synthesis. Proteins: Struct., Funct., Genet. 1995, 21, 167−195. (48) Onuchic, J. N.; Socci, N. D.; LutheySchulten, Z.; Wolynes, P. G. Protein Folding Funnels: The Nature of the Transition State Ensemble. Folding Des. 1996, 1, 441−450. (49) Onuchic, J. N.; LutheySchulten, Z.; Wolynes, P. G. Theory of Protein Folding: The Energy Landscape Perspective. Annu. Rev. Phys. Chem. 1997, 48, 545−600. (50) Koga, N.; Takada, S. Roles of Native Topology and ChainLength Scaling in Protein Folding: A Simulation Study with a Go-Like Model. J. Mol. Biol. 2001, 313, 171−180.

(51) Miyazawa, S.; Jernigan, R. L. Residue-Residue Potentials with a Favorable Contact Pair Term and an Unfavorable High Packing Density Term, for Simulation and Threading. J. Mol. Biol. 1996, 256, 623−644. (52) Chu, J. W.; Voth, G. A. Allostery of Actin Filaments: Molecular Dynamics Simulations and Coarse-Grained Analysis. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 13111−13116. (53) Lyman, E.; Pfaendtner, J.; Voth, G. A. Systematic Multiscale Parameterization of Heterogeneous Elastic Network Models of Proteins. Biophys. J. 2008, 95, 4183−4192. (54) Zhang, Z. Y.; Lu, L. Y.; Noid, W. G.; Krishna, V.; Pfaendtner, J.; Voth, G. A. A Systematic Methodology for Defining Coarse-Grained Sites in Large Biomolecules. Biophys. J. 2008, 95, 5073−5083. (55) 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. (56) 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. (57) Reinstein, J.; Schlichting, I.; Wittinghofer, A. Structurally and Catalytically Important Residues in the Phosphate Binding Loop of Adenylate Kinase of Escherichia-Coli. Biochemistry 1990, 29, 7451− 7459. (58) 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. (59) 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. (60) 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.; et al. Intrinsic Motions Along an Enzymatic Reaction Trajectory. Nature 2007, 450, 838−844. (61) 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. (62) 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. (63) Daily, M. D.; Phillips, G. N.; Cui, Q. Many Local Motions Cooperate to Produce the Adenylate Kinase Conformational Transition. J. Mol. Biol. 2010, 400, 618−631. (64) Wang, Y.; Gan, L. F.; Wang, E. K.; Wang, J. Exploring the Dynamic Functional Landscape of Adenylate Kinase Modulated by Substrates. J. Chem. Theory Comput. 2013, 9, 84−95. (65) Das, A.; Gur, M.; Cheng, M. H.; Jo, S.; Bahar, I.; Roux, B. Exploring the Conformational Transitions of Biomolecular Systems Using a Simple Two-State Anisotropic Network Model. PLoS Comput. Biol. 2014, 10, e1003521. (66) 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. (67) Gotz, A. W.; Williamson, M. J.; Xu, D.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simulations with Amber on Gpus. 1. Generalized Born. J. Chem. Theory Comput. 2012, 8, 1542−1555. (68) Salomon-Ferrer, R.; Gotz, A. W.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simulations with Amber on Gpus. 2. Explicit Solvent Particle Mesh Ewald. J. Chem. Theory Comput. 2013, 9, 3878−3888. (69) Zuckerman, D. M. Simulation of an Ensemble of Conformational Transitions in a United-Residue Model of Calmodulin. J. Phys. Chem. B 2004, 108, 5127−5137. (70) Best, R. B.; Chen, Y.-G.; Hummer, G. Slow Protein Conformational Dynamics from Multiple Experimental Structures: I

DOI: 10.1021/acs.jpcb.5b07668 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B The Helix/Sheet Transition of Arc Repressor. Structure 2005, 13, 1755−1763. (71) Okazaki, K. i.; Koga, N.; Takada, S.; Onuchic, J. N.; Wolynes, P. G. Multiple-Basin Energy Landscapes for Large-Amplitude Conformational Motions of Proteins: Structure-Based Molecular Dynamics Simulations. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 11844−11849. (72) Zheng, W. J.; Brooks, B. R.; Hummer, G. Protein Conformational Transitions Explored by Mixed Elastic Network Models. Proteins: Struct., Funct., Genet. 2007, 69, 43−57. (73) Takagi, F.; Kikuchi, M. Structural Change and Nucleotide Dissociation of Myosin Motor Domain: Dual Go̅ Model Simulation. Biophys. J. 2007, 93, 3820−3827. (74) Okazaki, K. i.; Takada, S. Dynamic Energy Landscape View of Coupled Binding and Protein Conformational Change: Induced-Fit Versus Population-Shift Mechanisms. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 11182−11187. (75) Tama, F.; Gadea, F. X.; Marques, O.; Sanejouand, Y. H. Building-Block Approach for Determining Low-Frequency Normal Modes of Macromolecules. Proteins: Struct., Funct., Genet. 2000, 41, 1− 7. (76) Lu, M.; Ma, J. The Role of Shape in Determining Molecular Motions. Biophys. J. 2005, 89, 2395−2401. (77) Ponder, J. W.; Case, D. A. Force Fields for Protein Simulations. Adv. Protein Chem. 2003, 66, 27−85. (78) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Struct., Funct., Genet. 2006, 65, 712−725. (79) 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. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2015, 5, 310−323. (80) 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. (81) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules 0.1. The Method. J. Comput. Chem. 1992, 13, 1011−1021. (82) Li, W. F.; Terakawa, T.; Wang, W.; Takada, S. Energy Landscape and Multiroute Folding of Topologically Complex Proteins Adenylate Kinase and 2ouf-Knot. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 17789−17794. (83) Le Grand, S.; Götz, A. W.; Walker, R. C. Spfp: Speed without Compromisea Mixed Precision Model for Gpu Accelerated Molecular Dynamics Simulations. Comput. Phys. Commun. 2013, 184, 374−380. (84) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald - an N.Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (85) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (86) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. Molecular-Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (87) Beckstein, O.; Denning, E. J.; Perilla, J. R.; Woolf, T. B. Zipping and Unzipping of Adenylate Kinase: Atomistic Insights into the Ensemble of Open↔Closed Transitions. J. Mol. Biol. 2009, 394, 160− 176. (88) Gur, M.; Madura, J. D.; Bahar, I. Global Transitions of Proteins Explored by a Multiscale Hybrid Methodology: Application to Adenylate Kinase. Biophys. J. 2013, 105, 1643−1652. (89) Sugita, Y.; Okamoto, Y. Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 1999, 314, 141−151.

(90) Mitsutake, A.; Sugita, Y.; Okamoto, Y. Generalized-Ensemble Algorithms for Molecular Simulations of Biopolymers. Biopolymers 2001, 60, 96−123.

J

DOI: 10.1021/acs.jpcb.5b07668 J. Phys. Chem. B XXXX, XXX, XXX−XXX