Free-Energy Landscape of Protein–Ligand Interactions Coupled with

Jan 4, 2017 - (B) Ligand-binding pocket of the GlnBP–glutamine complex. The glutamine ligand is shown in green. The residues in the large domains ar...
0 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCB

Free-Energy Landscape of Protein−Ligand Interactions Coupled with Protein Structural Changes Kei Moritsugu,*,† Tohru Terada,‡ and Akinori Kidera† †

Graduate School of Medical Life Science, Yokohama City University, 1-7-29 Suehirocho, Tsurumi, Yokohama 230-0045, Japan Graduate School of Agricultural and Life Sciences, The University of Tokyo, 1-1-1 Yayoi, Bunkyo, Tokyo 113-8657, Japan

‡

S Supporting Information *

ABSTRACT: Protein−ligand interactions are frequently coupled with protein structural changes. Focusing on the coupling, we present the free-energy surface (FES) of the ligand-binding process for glutamine-binding protein (GlnBP) and its ligand, glutamine, in which glutamine binding accompanies large-scale domain closure. All-atom simulations were performed in explicit solvents by multiscale enhanced sampling (MSES), which adopts a multicopy and multiscale scheme to achieve enhanced sampling of systems with a large number of degrees of freedom. The structural ensemble derived from the MSES simulation yielded the FES of the coupling, described in terms of both the ligand’s and protein’s degrees of freedom at atomic resolution, and revealed the tight coupling between the two degrees of freedom. The derived FES led to the determination of definite structural states, which suggested the dominant pathways of glutamine binding to GlnBP: first, glutamine migrates via diffusion to form a dominant encounter complex with Arg75 on the large domain of GlnBP, through strong polar interactions. Subsequently, the closing motion of GlnBP occurs to form ligand interactions with the small domain, finally completing the native-specific complex structure. The formation of hydrogen bonds between glutamine and the small domain is considered to be a rate-limiting step, inducing desolvation of the protein−ligand interface to form the specific native complex. The key interactions to attain high specificity for glutamine, the “door keeper” existing between the two domains (Asp10−Lys115) and the “hydrophobic sandwich” formed between the ligand glutamine and Phe13/Phe50, have been successfully mapped on the pathway derived from the FES.

■

INTRODUCTION The complex biological networks of biochemical reactions and transduction processes are frequently regulated by interactions of ligands with receptor proteins, which are accompanied by various levels of structural changes of the proteins. Computational studies based on all-atom molecular dynamics (MD) simulations have made significant contributions toward understanding such dynamical processes of ligand binding, which include long-time simulations,1−3 the metadynamics method,4,5 and analyses based on the Markov state model.6−10 These results provided various pictures of ligand-binding processes, covering diffusion on the protein surface, the formation of a transient “encounter complex”, and recognition at the ligandbinding site. In studies of ligand-binding processes, the coupling between protein structural change and ligand binding is one of the most important subjects, as has been discussed in terms of the two classical models, induced fit11 and conformational selection.12 However, the complexity lies in description of the coupling at the atomistic-interaction level: During the binding process, how do the dynamical structures of the protein modulate and regulate the ligand-binding modes? Conversely, how do the ligand interactions influence the ensemble of protein structures? The difficulty in obtaining such descriptions © XXXX American Chemical Society

with sufficient statistics is due partly to the time scale separation among ligand diffusion, ligand−protein interaction, and protein structural change, which causes complications in capturing the different types of dynamics in a simulation process. To circumvent this difficulty in time scale separation, we adopted a method sampling configurational space, the multiscale enhanced sampling (MSES) method,13−17 to calculate the free-energy landscape in both the protein’s and ligand’s degrees of freedom, focusing on the atomistic protein−ligand association after ligand diffusion or the process starting from the encounter complex and ending at the native bound state. MSES is used to enhance configurational sampling by adopting a multiscale scheme, which allows coupling of an all-atom model (MM) with the accelerated dynamics of the coarsegrained (CG) degrees of freedom defined by prior knowledge or experimental data, in combination with the Hamiltonian replica exchange18 that eliminates the bias via coupling to the CG model. The difficulty in application to large protein systems, or the scalability problem in enhanced sampling, was Received: November 20, 2016 Revised: December 27, 2016 Published: January 4, 2017 A

DOI: 10.1021/acs.jpcb.6b11696 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

has been well studied both experimentally and theoretically as a model protein−ligand complex.19,21−23 The crystal structure of the glutamine-bound GlnBP revealed that the high ligand specificity is attained not only by polar interactions with the glutamine ligand in the binding pocket (Table 1) but also by the “door keeper”, the interdomain

solved by limiting the conformational sampling of MM to the “essential subspace”, which is suitably defined with small dimensions of the CG model (see Methods). Our previous studies on the folding dynamics of chignolin,13 the ordering transition of an intrinsically disordered protein (sortase),14 and the protein−protein interaction of the barnase−barstar complex16 have exhibited a wide range of applicabilities of enhanced samplings to large protein systems on atomistic resolution and in explicit solvent. For the simulation system of ligand-binding dynamics coupled with protein structural change, we chose glutaminebinding protein (GlnBP) and its substrate, glutamine. GlnBP is a periplasmic binding protein in Gram-negative bacteria, and it is responsible for glutamine permeation across the cell membrane. Both the ligand-bound19 and unbound20 structures were solved by X-ray crystallography and showed that a ligandinduced structural change occurs in the form of bending of the hinge (residues 85−89 and 181−185) connecting the two globular domains, “large domain” (residues 1−84 and 186− 226) and “small domain” (residues 90−180) (Figure 1A).19 The motion occurring during ligand binding is a typical domain closure, leading to highly specific recognition. This dynamic recognition is common to other periplasmic binding proteins, like histidine-binding protein and lysine, arginine, ornithine (LAO) protein, and the smallest system among them, GlnBP,

Table 1. Probability of Hydrogen-Bond Formation, pHB, at the Protein−Ligand Interface during the Simulations pHB,MSESc #

protein

ligand

site

pHB,MM

1 2 3 4 5 6 7 8 9 10

Asp10Oδ Ala67O Thr70N Arg75Nη Lys115Nε Gly119N Asp157Oδ Gly68O Thr70Oγ2 His156Nε2

GlnNε2 GlnNε2 GlnO GlnO GlnOε1 GlnO GlnN GlnN GlnN GlnOε1

L1,2

0.99 1.00 0.99 1.00 0.99 1.00 1.00 0.51 0.67 0.76

a

L3,4 S5,6,7

b

S1

S2

S3

0.99 1.00 0.99 1.00 0.99 1.00 0.85 0.83 0.86 0.81

0.63 0.65 0.66 0.83 0.37 0.53 0.33 0.41 0.47 0.04

0.48 0.57 0.70 0.90 0.07 0.14 0.03 0.35 0.31 0.01

a

A list of the hydrogen bonds was made according to those formed in the bound/closed crystal structure.19 Sites 1−7 are those stably maintained in the 50 ns MM simulation starting from the bound/ closed crystal structure. These hydrogen bonds were used for the descriptions of substates, S1, S2, and S3. bpHB,MM represents the probability of occurrence of hydrogen bonds in the MM simulation. c pHB,MSES represents the probability in the unbiased MSES ensemble. S1 if the door keeper between Asp10 and Lys115 is formed. Otherwise, S2 if RMSDprot < 6 Å. Otherwise, S3.

hydrogen bonds between Asp10 (in the large domain) and Lys115 (in the small domain), and “hydrophobic sandwich”, a hydrophobic pocket formed by Phe13 and Phe50 (Figures 1B and S1).19 Similar amino acids, including glutamate, asparagine, and aspartate, showed a low binding affinity toward GlnBP.24 To scrutinize the dynamical process that establishes such subtle interactions, detailed atomic information on the large configurational space is crucially important. Here, the structural ensemble for the formation of the GlnBP−glutamine complex, ranging from the encounter complex to the natively bound state, was fully calculated by the MSES simulation and analyzed with the protein−ligand interactions and magnitude of protein structural changes as the reaction coordinates. These studies illustrated the energy landscape at atomic resolution and included the surrounding solvent molecules, which clarified the mechanism determining the coupling behavior between ligand interaction and structural change of the protein.

■

METHODS MSES Simulation. Here, we explain the MSES method briefly. The detailed formulas are described in the literature.13,15,17 The simulation system comprises an MM of protein and ligand molecules with surrounding solvents (rMM; N degrees of freedom) and the associated CG model (rCG; M degrees of freedom). The potential energy of the multiscale system is given by

Figure 1. (A) Glutamine-bound19 (green) and unbound20 (gray) crystal structures of GlnBP, shown after superimposition of the large domain. (B) Ligand-binding pocket of the GlnBP−glutamine complex. The glutamine ligand is shown in green. The residues in the large domains are shown in red, those in the small domains are shown in blue, and the two phenylalanines (Phe13 and Phe50) related to the hydrophobic sandwich are shown in orange.

V = VMM(rMM) + VCG(rCG) + VMMCG(rMM, rCG, kMMCG) (1)

with B

DOI: 10.1021/acs.jpcb.6b11696 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B VMMCG = kMMCG[χMM (rMM) − χCG (rCG)]2

replicas.18 The exchange probability between replicas m and n, so as to satisfy the detailed balance conditions, is given at a high-temperature limit of adiabatic separation (eqs 3 and 7) by17

(2)

where VMM and VCG are the MM and CG potential energy functions, respectively. The term VMMCG is the coupling (harmonic constraint) for K collective variables, χCG(rCG), which is defined by CG coordinates. χMM(rMM) is also a Kdimensional vector, which is a projection of rMM onto the associated K-dimensional space. The coupling constant, kMMCG, is used to tune the acceleration of the MM dynamics via coupling with CG. To derive the free-energy surface (FES) solely from VMM, without bias from coupling VMMCG, we calculated the joint distribution (β is the inverse temperature of the system) and extrapolated it to kMMCG = 0 in the process of the MSES simulation. The joint distribution was derived by a simulation protocol of Gibbs sampling; that is, by sampling rMM, rCG, and kMMCG separately in an iterative manner with the conditional probabilities.25,26 For efficient samplings of rMM and rCG, the approximation of adiabatic separation was introduced by making CG much heavier than MM, to move much slower than MM, while making the CG temperature, 1/β′, much higher than that of MM, that is, β′ ≪ β.27−33 Under the approximation, the conditional probabilities for rMM and rCG are written respectively as28,32

pmn = min(1, exp(Δmn))

with n m m n ⎡ ρ(kMMCG |r mMM, r CG |r nMM, r CG )ρ(kMMCG )⎤ Δmn = ln⎢ ⎥ m m n n |r mMM, r CG |r nMM, r CG )ρ(kMMCG )⎦ ⎣ ρ(kMMCG m n m = β(kMMCG − kMMCG )([χMM (r mMM) − χCG (r CG )]2 n − [χMM (r nMM) − χCG (r CG )]2 ) m n n m ⎡ Z(r CG, kMMCG)Z(r CG, kMMCG )⎤ − ln⎢ ⎥ m m n n ⎣ Z(r CG, kMMCG)Z(r CG, kMMCG) ⎦ (9)

It is noted that the second term on the right side of eq 9 is not easy to evaluate. To solve this problem, we utilize the condition in eq 7 that CG does not encounter the influence from the MM system, allowing the replicated systems of the MCMC simulation to consist of many MMs driven by a single copy of CG, that is, rmCG = rnCG. This system setup allows the numerator and denominator of the second term on the right side of eq 9 to be canceled out, and we finally have, instead of eq 9

ρ(rMM|rCG, kMMCG)

m n m Δmn = β(kMMCG − kMMCG )([χMM (r mMM) − χCG (r CG )]2

= exp( −β[VMM(rMM) + VMMCG(rMM, rCG, kMMCG)]) /Z(rCG, kMMCG)

n − [χMM (r nMM) − χCG (r CG )]2 )

(3)

(4)

with Z(rCG, kMMCG) ≡

∫ drMM

exp( −β[VMM(rMM) + VMMCG(rMM, rCG, kMMCG)])

Z(kMMCG) ≡

(5)

∫ exp[−β′VCG(rCG)]Z(rCG, kMMCG)β′/β drCG (6)

In eq 4, rMM does not appear explicitly in the argument because the slow motion of CG allows preaveraging of MM motions, and the influence of MM is counted implicitly in the potential of mean force, Z(rCG, kMMCG)β′/β.28,32 Here, we further take the high-temperature limit in eq 4, or β′/β → 0,17 that is ρ(rCG|kMMCG) = exp[−β′VCG(rCG)]/Z(kMMCG)

(10)

In summary, the MSES simulation process consists of iterated processes, MD simulations for rMM (eq 3) and rCG (eq 7), followed by MCMC for kMMCG using eqs 8 and 10. The smallness of K, or the dimensions of χMM and χCG that are determined by the CG degrees of freedom, guarantees a small value of Δmn or a large exchange probability pmn, which allows applications to large protein systems. Potential Energy Functions and Kinetic Parameters for MSES. Both the MM and CG energy functions were determined as follows. AMBER ff99SBildn was used for the allatom potential energy, VMM.34 To focus on the coupling between ligand binding and protein structural change, GlnBP and glutamine were constrained so that they were not too far from each other, using the two constraint forces between the Cα atoms of residues 68/118 at the GlnBP large/small domains and the Cα atom of the glutamine ligand (force constant = 100 kcal/mol/Å2 when d > 20 Å). The CG model consists of the Cα atoms of GlnBP and glutamine; that is, M = (226 + 1) atoms × 3. The CG potential energy, VCG, was set as the sum of two terms for representing both the intraprotein and protein− ligand interactions. For the intraprotein interactions, the mixed elastic network model was adopted,35 by embedding the bound (PDB: 1WDN)19 and unbound (PDB: 1GGG)20 structures as the minima of the two potential wells. The force constant and cutoff length were set to 1.8 kcal/mol/Å2 and 12 Å in the elastic network model. The mixing temperature was 2 × 106 K, and the energy difference, ΔE, of the two basins was 55 kcal/mol larger in the unbound state than that in the bound state. For protein−ligand interactions, the 6-12 Lennard-Jones potential was applied to 84 pairs of Cα atoms between GlnBP and glutamine, together with a soft harmonic boundary to avoid an overly large separation between the protein and ligand, as in our previous study of protein−protein interactions.16 The 84 Cα atoms of GlnBP were selected as those located within 12 Å

ρ(rCG|kMMCG) = exp[−β′VCG(rCG)]Z(rCG, kMMCG)β ′ / β /Z(kMMCG)

(8)

(7)

It is noted that the denominators in eqs 3 and 7 can be ignored during the MD simulation because they are simply taken as the normalization constants. Equations 3 and 7 produce the largest driving force for MM because CG can move freely without experiencing the counterforce from the potential of mean force of MM.17 The temperature and mass of the CG system were set to satisfy the conditions of adiabatic separation and β′/β → 0, according to ref 17. Sampling in terms of kMMCG was then performed by Markov chain Monte Carlo (MCMC), for discretized values of kMMCG. Here, the Hamiltonian replica exchange was used, where many replicated systems with various kMMCG values, ranging from a large value to zero, are exchanged between neighboring C

DOI: 10.1021/acs.jpcb.6b11696 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

asparagine, glutamate and aspartate were performed, starting at the bound structure for 50 ns and under the same simulation conditions as those described above.

from the Cα atom of the ligand in the crystal structure of the bound state. The minimum depth and distance for the LJ potential were set as 0.04 kcal/mol and the distance in the crystal structure of the bound state, respectively. The soft boundary potential starts at 12 Å apart from the minimum of the LJ potential, with a harmonic force constant of 50 kcal/ mol/Å2. Under the conditions of adiabatic separation, the mass of the CG atom must be sufficiently large, and it was set to 2 × 104. For the coupling potential, VMMCG, 113 pairs of Cα atoms between the large and small domains of GlnBP and 84 pairs of Cα atoms between GlnBP and glutamine were used as the variables of χMM and χCG (i.e., K = 197 and therefore K < M ≪ N ∼ 105). These pairs were chosen according to the condition that their distances were less than 12 Å in the crystal structure of the bound state. Computations. The starting structure for the simulations was taken from the crystal structure of the bound state, and the missing N- and C-terminal residues were modeled by MODELLER.36 Note that an oxygen atom in the main chain of ligand glutamine is lacking in the associated PDB file and thus it was also modeled (the modeled structure is referred to as the “crystal structure” in this study). The MSES simulations were carried out using the class library for multicopy and multiscale MD simulations, μ2lib.37 The MM was constructed in a rectangular box, with a margin of 12 Å to the simulation box boundary, containing ∼12 000 TIP3P water molecules38 together with two sodium ions to neutralize the system, which yields a total of 39 425 atoms. MM simulations were performed with a time step of 2 fs using the SHAKE algorithm,39 constraining the bonds involving hydrogen atoms; under constant pressure and temperature (NPT) conditions of P = 1 atm and T = 300 K, maintained using Berendsen’s thermostat and a barostat;40 and at a relaxation time of 1 ps. Electrostatic interactions were calculated with the particle-mesh Ewald method.41 The CG simulation was performed using Berendsen’s thermostat under the constant temperature (NVT) condition at T = 3 × 103 K and a mass of CG atoms of 2 × 104. It is shown in Figure S2 that the CG motion was sufficiently slow so that the condition of adiabatic separation between MM and CG was satisfied in the present simulation system, whereas a significantly smaller CG mass, 2 × 102, accelerated the CG motions to break the approximation. The MM simulation was performed under the restraints from CG (eq 3). The slowness of the CG particles allowed the CG coordinates in eq 3 to be updated every 20 ps.17 Here, MCMC was performed every 20 ps, using 12 replicas, with kMMCG = 0, 0.00032, 0.0014, 0.0038, 0.0078, 0.0136, 0.0215, 0.0318, 0.0447, 0.0603, 0.0789, and 0.1007 kcal/mol/Å2. These values were obtained by the relation log(VMMCG) = −0.23 × log(kMMCG) + 5.2, which was evaluated by test calculations with a single replica for various kMMCG values, while referring to eq 8 to make the exchange rates of all neighboring replica pairs equal to each other. Five MSES simulation runs were performed for 50 ns each, amounting to 12 × 50 ns × 5 runs = 3 μs in total. For the purpose of having sufficient samples around the bound state, the simulation trajectories were omitted in the analysis of the unbiased MSES ensembles, when they stayed in the unbound state (RMSDprot > 12 Å and dlig > 15 Å) for over 1 ns, according to the procedure of ref 42. Hence, the trajectory lengths used for the analysis were 50, 50, 44, 44, 40 ns for the five runs, respectively (see Figure S3). For comparison, conventional equilibrium MD runs of the wild type as well as the mutants and GlnBP complexed with

■

RESULTS AND DISCUSSION MSES Simulation and Free-Energy Landscape. An MSES simulation of five 50 ns runs for the GlnBP−glutamine complex in explicit solvent was carried out to fully sample the atomistic configurations during the ligand-binding process coupled with the structural change in GlnBP. Twelve replicas for the Hamiltonian replica exchange were found to be sufficient for simulating the solvated system containing ∼40 000 atoms, owing to the highly scalable MSES.13−17 The distributions of the MM/CG coupling energies, VMMCG (eq 2 in Methods), of the 12 replicas significantly overlapped for neighboring replicas (Figure S3A), guaranteeing a high exchange probability; the average acceptance ratio of all exchanges for the five 50 ns runs was 0.24, indicating successful Hamiltonian exchange simulation. The CG model, VCG, used in this study consists of the mixed elastic network model,25 for the intraprotein interactions generating transitions between the closed and open forms, and the Lennard-Jones-type potential, for the protein−ligand interactions (see Methods for details). The ligand molecule is also represented by a single Cα atom. The parameters for VCG were designed to efficiently accelerate the structural transitions. The values of RMSDprot (the Cα root-mean-square displacement from the bound/closed crystal structure) measuring the protein structural change (Figure S3B) indicated that the distribution obtained from the CG simulation (the CG potential only, without coupling to the MM model) has only a small probability of occurrence in the bound/closed state, but the maximum exists in between the two terminal structures, thus accelerating frequent structure transitions. The strong coupling with the CG model through a large value of kMMCG in eq 1 can drive the MM to go over a large configurational space, reaching the ligand-unbound state. The FES was calculated from the unbiased MSES ensemble (VMM in eq 1), which was obtained by extrapolating kMMCG to zero in eq 1 by Hamiltonian exchange (eqs 8 and 10). In Figure 2A, the FES is represented in a two-dimensional (2D) space of RMSDprot and dlig [= dGly68−Gln (Cα of Gly68 in the large domain−glutamine ligand) + dThr118−Gln (Cα of Thr118 in the small domain−glutamine ligand)], representing the protein’s and ligand’s reaction coordinates, respectively, where dlig quantifies how close the ligand is to the bound/closed crystal structure. Note that the convergence of the simulation was confirmed in the comparison of the distributions for both RMSDprot and dlig, calculated from the trajectories of various lengths (Figure S4). These differences for 0−40 ns and for the full length, 0−50 ns, that is, |p (0−40 ns) − p (0−50 ns)|/p (0−50 ns), were found to be sufficiently small, with the averages being 0.057 and 0.046 for RMSDprot and dlig, respectively. Despite the narrow distribution in the MM simulation (the standard MM simulation with no use of the MSES scheme), starting from the bound/closed crystal structure, the distribution obtained by the MSES simulation was largely extended through coupling with the CG model. It is clarified that the unbiased MM in the MSES simulation (kMMCG = 0) does not follow the CG model completely but follows the MM force field to map the correct FES. The large variations in RMSDprot primarily represent the rigid-body hinge-bending motions of the two domains: the intradomain Cα RMSD values D

DOI: 10.1021/acs.jpcb.6b11696 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

process is quantified at the atomic interaction level by Q, a fraction of the occurrence of the 61 native contacts between GlnBP and the glutamine ligand observed in the bound/closed crystal structure (0 ≤ Q ≤ 1; the native contacts were defined by a distance less than 4 Å between nonhydrogen atoms, with more than 70% probability of occurrence in the present 50 ns MM simulation), as used in the protein folding and our previous protein−protein interaction studies.16 Figure 3A illustrates the FES on the 2D space of RMSDprot and Q, showing a funnel-like free-energy landscape going down from the unbound/open state to the native basin of the bound/ closed state (RMSDprot < 4 Å and Q > 0.8). However, even in this simplified 2D diagram, the FES does not exhibit a simple downhill landscape but a rugged surface, due to a bottleneck between the unbound/open state and the native bound/closed basin at RMSDprot = 6 or 7 Å and Q ∼0.5. This is characterized by off-pathway structures around RMSDprot < 5 Å and Q < 0.3, which are separated from the native basin by high free-energy barriers (see the next section for details). To elucidate the structural basis of the pathway linking the unbound/open state to the bound/closed state, as well as the origin of the associated bottleneck, on the basis of atomic interactions, we chose the pattern of hydrogen bonds between GlnBP and the glutamine ligand as a descriptor for the classification of the structural ensemble. Ten hydrogen bonds were found in the bound/closed crystal structure,19 and we calculated their probabilities of occurrence in a 50 ns MM simulation starting from the bound/closed crystal structure (Table 1). Among them, seven hydrogen bonds that were stably maintained in the MM simulation were chosen for the identification of the glutamine-bound state. In addition, as a hallmark of the native basin, we found a hydrogen bond between the two domains, formed between Asp10 (in the large domain) and Lys115 (in the small domain). In the study of the crystal structure, this hydrogen bond was named door keeper, due to the fact that it tethers the two domains and forms a front door that shields the binding pocket from the bulk solvent region (see Figures 1B and S1A).19 Prior to the classification based on atomic interactions, we first classify the structural ensemble roughly into the three distinct substates, S1, S2, and S3, in terms of the information on the GlnBP structure, that is, RMSDprot and the hallmark of the native basin, whether or not the door keeper is formed. State S1 describes the native basin, defined by the condition that the door keeper was formed. It was confirmed that state S1 overlapped with the native basin in the RMSDprot−Q space and in the pattern of the seven hydrogen bonds (Figure 3 and Table 1). State S2 was defined to satisfy the condition that the door keeper was not formed, but GlnBP still maintained the structure near the native region or RMSDprot < 6 Å. Despite the native-like protein structure, glutamine recognition in S2 is much weaker than that in S1, particularly at the small domain (see Table 1). This implies the importance of the door keeper in stabilizing the interactions between GlpBP and the ligand, as Asp10 and Lys115 bind the ligand only after the door keeper is formed. State S3 consists of structures with RMSDprot > 6 Å and with almost no probability of glutamine binding to the small domain. The pattern of the seven hydrogen bonds with the glutamine ligand (Table 1) further divided states S2 and S3 into five substates, S2LS, S2L, and S2U and S3L and S3U, respectively. The letter “L” designates the substate in which at least one hydrogen bond is formed at both sites L1,2 and L3,4, described in

Figure 2. (A) 2D FES plotted on RMSDprot vs. dlig plots. RMSDprot is the Cα root-mean-square displacement of the small domain from the bound/closed state after superimposing the large domain, and dlig is the sum of Cα distances, dGly68−Gln (Gly68 of the large domain− glutamine ligand) and dThr118‑Gln (Thr118 of the small domain-the ligand). (B) 2D FES plotted on dGly68−Gln vs. dThr118−Gln plots. The three FESs are given for the equilibrium all-atom MD simulation starting from the bound/closed structure (upper, MM), the CG simulation (middle, CG), and the unbiased ensemble derived from the MSES simulations (lower, MSES).

from the bound/closed crystal structure were less than 1.6 and 1.5 Å for the large and small domains, respectively, in all replicas. The positive correlation in the MSES distribution in Figure 2A suggests coupling between domain closure and ligand binding during the process of bound/closed structure formation. The FES is also illustrated by focusing on the binding of the ligand to the two domains in terms of the components of dlig, dGly68−Gln vs. dThr118−Gln (Figure 2B: their time courses are given in Figures S3C and S3D, with comparison to those of the MM simulation). The FES of the CG simulation differs greatly from that of the unbiased MSES ensemble because the CG model lacks electrostatic interactions. Ligand binding is initially driven by electrostatic interactions, particularly at the early stage of the binding process (see below). The large probability of occurrence in the region of small dGly68−Gln values suggests the existence of preferential electrostatic interactions at the large domain. Substates of the Ligand-Binding Process. Let us delve into the atomic details of the ligand-interaction processes observed in the unbiased MSES ensemble. Here, instead of dlig using only the coordinates of the Cα atoms, the ligand-binding E

DOI: 10.1021/acs.jpcb.6b11696 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 3. (A) 2D free-energy surfaces along RMSDprot and Q for the unbiased ensemble derived from the MSES simulations. (B) Mapping the substates on a pathway from the unbound/open state (S3U) to the bound/closed state (S1), with representative structures of the small domain (colored ribbons) and Cα atoms of the glutamine ligand (magenta dots). These were drawn after superimposing the large domain. Each portion of the 2D FES along RMSDprot and Q is also shown for each substate.

Table 2. Probabilities and Quantities of Protein Structure and Ligand Binding for the Six Substates substate

S1

S2LS

S2L

S2U

S3L

S3U

pa Qb RMSDprot [Å]b RMSDlig [Å]b,c nHB for L1,2b,d nHB for L3,4b,d nHB for S5,6,7b,d ddoor keeper [Å]b,e nCb,f nC,Pheb,g nW,protb,h nW,ligb,i

1 0.85 ± 0.09 3.1 ± 0.83 0.72 ± 0.30 2.0 ± 0.058 2.0 ± 0.087 2.8 ± 0.36 2.8 ± 0.13 71 ± 6.0 20 ± 4.8 16 ± 2.2 1.9 ± 0.76

0.15 0.62 ± 0.17 4.3 ± 0.92 1.4 ± 0.59 1.9 ± 0.33 1.9 ± 0.26 2.2 ± 0.74 4.7 ± 0.85 57 ± 9.1 14 ± 6.5 21 ± 3.8 4.2 ± 1.9

0.053 0.35 ± 0.10 5.2 ± 1.5 2.3 ± 1.0 1.9 ± 0.25 1.9 ± 0.32 0 7.3 ± 0.97 39 ± 8.5 8.4 ± 6.7 34 ± 5.0 8.8 ± 2.1

0.10 0.13 ± 0.14 4.5 ± 1.8 5.9 ± 2.5 0.70 ± 0.80 0.046 ± 0.29 0.39 ± 0.70 6.4 ± 1.4 29 ± 9.2 5.7 ± 4.3 30 ± 3.7 10 ± 2.2

0.51 0.36 ± 0.12 10.4 ± 3.5 2.5 ± 1.3 1.8 ± 0.37 1.8 ± 0.41 0 10 ± 3.3 40 ± 9.2 8.7 ± 7.0 33 ± 5.1 9.0 ± 2.2

0.36 0.18 ± 0.13 10.2 ± 3.9 6.5 ± 2.8 1.3 ± 0.81 0.028 ± 0.21 0.15 ± 0.75 10 ± 3.2 28 ± 10 5.0 ± 4.6 35 ± 4.7 12 ± 3.1

a

Probability of occurrence of each substate, p, relative to p(S1). Note that the p value of S3U is underestimated due to the simulation analysis used in this study, in which completely dissociated states were ignored (see Methods for details). bAverage and standard deviation of various structural quantities within the substate. cRMSDlig is the RMSD value for the nonhydrogen atoms of the ligand from the bound/closed position, after superimposition on large domains. dnHB is the number of hydrogen bonds formed for L1,2, L3,4, or S5,6,7. The maxima of these values are 2, 2, and 3, respectively. eddoor keeper is the distance between Oδ of Asp10 and Nζ of Lys115, constituting the door keeper. fnC is the number of contacts between GlnBP and the ligand, including non-native ones. gnC,Phe is the number of contacts between Phe13 and Phe50, constituting the hydrophobic sandwich and the ligand. hnW,prot is the number of hydrated waters within 4 Å from the nonhydrogen atoms of the binding pocket of GlnBP. The pocket is defined by the nonhydrogen atoms of GlnBP interacting with the ligand in the bound/closed state. inW,lig is the number of hydrated waters within 4 Å from the nonhydrogen atoms of the ligand.

Table 1. Thus, at least two hydrogen bonds are formed between the large domain and ligand atoms. The letter “S” represents

binding with at least one hydrogen bond to site S5,6,7, defined in Table 1. “LS” indicates that events L and S occur concurrently, F

DOI: 10.1021/acs.jpcb.6b11696 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B and “U” indicates that neither L nor S occurs, although U can have ligand contacts other than those in L or S (i.e., Q ≥ 0). As the occurrence of S2S and S3S was found to be negligible, these substates were classified to S2U and S3U, respectively. In summary, six substates, S1, S2LS, S2L, S2U, S3L, and S3U, were identified during the ligand-binding process (see Figure 3). The associated quantities characterizing the structures of the six substates are summarized in Table 2. The distributions of these quantities are also given in Figure S5, indicating the definitions of the substates that are clearly distinct from each other in terms of ligand interaction and protein structure. Pathway Connecting the Substates. The mutual disposition of the distributions of the substates (see Figure 3), as well as the various structural characteristics (see Table 2), particularly the number of hydrogen bonds with the ligand, ⟨nHB⟩, suggests the dominant pathway of the ligand-binding process to be S3U−S3L−S2L−S2LS−S1, whereas S2U is the off-pathway substate, exhibited in Figure 3B. At the start of the binding process, a glutamine migrates via diffusion to interact with a binding site of the large domain in the open form (S3U) (Arg75/Thr70 at site L1,2, which is closer to the hinge region than the other site, L3,4, with Asp10/Ala67, see Figure 1B), mostly by electrostatic interactions between the positive charge of the Arg75 side chain and the negative charge of the mainchain carboxyl group. These features in S3U are clearly seen in the value of ⟨nHB⟩ for L1,2, which is greater than 1 (see Table 2). However, we noted that S3U contains numerous misbound structures with large RMSDlig values, the RMSD value for the ligand from the native bound/closed structure (non-native interactions not included in Table 1, see below). At the next stage, there are two directions. The first is the on-pathway route to S3L, which allows the glutamine to enhance the interactions with sites L1,2 and L3,4 further (⟨nHB⟩ in Table 2 close to 2) and stabilize the ligand’s position at the large domain (decrease in ⟨ΔRMSDlig2⟩1/2 from 6.5 Å (S3U) to 2.5 Å (S3L)). At S3L, the stable interactions with the large domain enable two phenylalanine residues in the large domain, Phe13 and Phe50, to start packing the glutamine ligand (increase in ⟨nC,Phe⟩ from 5.0 (S3U) to 8.7 (S3L) in Table 2, where nC,Phe is the number of contacts between the two phenylalanine residues and the ligand). In the study of the crystal structure, this interaction was named hydrophobic sandwich19 (Figures 1B and S1B), contributing to ligand binding and stabilization of the small loop of residues 11−16, containing Phe13.43 Substate S3L resembles the encounter complex, “state 11”, found in the Markov state analysis of another periplasmic binding protein, LAO binding protein, in which the ligand interacts with two aromatic residues, Tyr14 and Phe52.7 The other is the offpathway route, which causes the closing motion of the two domains before the glutamine properly binds (S2U). S2U is characterized by almost closed domain structures (small RMSDprot, 4.5 Å) and weak interactions with the ligand (large RMSDlig, 5.9 Å and small Q, 0.13). After S3L, along the on-pathway route to S2LS, closing motion of the two domains occurs to form ligand interactions at site S5,6,7 (shown by the small RMSDprot of 4.3 Å and the large nHB for S5,6,7 of 2.24). Finally, the door keeper, or the salt bridge between Asp10 and Lys115, closes the binding cleft to complete the ligand-bound structure at S1 or the native complex structure. The two residues, Asp10 and Lys115, already start to approach each other in substate S2LS (ddoor keeper decreases from 7.3 Å (S2L) to 4.7 Å (S2LS), where ddoor keeper is the distance between

Asp10 and Lys115). Therefore, substate S2LS is on the pathway toward S1. The binding process is also viewed through the change in the fluctuations of the small domain and glutamine ligand (the standard deviations of RMSDprot and RMSDlig, respectively, in Table 2). The decrease in the fluctuations upon closing motions of GlnBP and glutamine binding demonstrates again the funnel-like landscape of the ligand-binding process, except for off-pathway substate S2U, which retains large fluctuations for both the small domain and ligand. In substates S2U and S3U, the ligand frequently binds in non-native modes, as listed in Table S1. These binding modes often lead to misbound structures, such as the example shown in Figure S6, in which the ligand binds in the opposite direction. In this structure, all four non-native binding modes described in Table S1 exist. However, the non-native modes are completely filtered out before the subsequent substate, S3L, in which only the binding mode satisfying the native interactions at both sites, L1,2 and L3,4, is allowed. Solvation and the Rate-Limiting Step. The ligandbinding process was further illustrated in terms of solvation (nW,prot and nW,lig in Table 2 and Figure S5, where nW,prot is the number of hydrated waters in the binding pocket of GlnBP and nW,lig is the number of hydrated waters bound to the ligand). The results clearly demonstrated that desolvation of the binding sites occurs upon the closing motion or upon formation of substate S2LS (decreases in nW,prot and nW,lig from 34 and 8.8 (S2L) to 21 and 4.2 (S2LS), respectively), together with subsequent formation of the door keeper, suggesting that the rate-limiting step of the ligand-binding process may exist between S2L and S2LS. This inference can be qualitatively verified by the following MM simulations. In the MM simulation starting from the native bound/closed state, or S1, the protein stayed at S1 (Figures 2 and S3). Although a certain population reached S2LS (the probability of occurrence was 0.06), the simulation never reached S2L within 50 ns. This indicates a possibility of an energy barrier existing between S2LS and S2L. In other words, even though the door keeper is transiently disrupted, the protein is able to maintain the closed structure. To further examine the possible energy barrier between S2LS and S2L, additional equilibrium MM simulations were performed, starting from S2LS (a 50 ns simulation) and S2L (two 50 ns simulations). The starting structures were chosen from the trajectory satisfying the condition that they are sufficiently close to the boundary of the two substates on the RMSDprot−Q space (see Figure 3). Figure 4A clearly shows that the simulation starting from S2LS goes back to S1, whereas the simulations starting from S2L reach S3. It is thus considered that the ratelimiting step during the ligand-binding process may exist between S2LS and S2L, where the GlnBP closing motion causes the desolvation of the GlnBP−glutamine interface as well as the association of the glutamine ligand with the small domain. Roles of Ligand-Binding Sites. To examine the roles of sites L1,2, L3,4, and S5,6,7 in the ligand-binding process, MM simulations of three mutants, R75A (L1,2), D10A (L3,4), and K115A (S5,6,7), eliminating the ligand-binding side chain at each site, were further performed, starting from the bound/closed crystal structure, for 50 ns. The simulation results indicated that R75A that is mutated on a residue at L1,2 lost the ligand interactions at all three sites, leading to domain opening (Figure 4B) and demonstrating the importance of the polar G

DOI: 10.1021/acs.jpcb.6b11696 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

5A). In the simulations, Asn did not remain stable at the binding sites of either the large or small domain and caused

Figure 4. (A) 2D profiles of the distribution on RMSDprot and Q for 50 ns equilibrium all-atom MD simulations, one starting from a structure at S2BL and the others from two different structures belonging to S2L. (B) RMSDprot and the minimum distances for hydrogen bonds in L1,2 (dHB,L1,2), L3,4 (dHB,L3,4), S5,6,7 (dHB,S5,6,7), and in #1 and 5 (related to door-keeper formation, dHB,1,5) are plotted as a function of time. Equilibrium all-atom MD simulations (50 ns) of the wild type (MM) and three mutants, D10A, R75, and K115A.

contact at Arg75 as the initial binding site occurring at S3U. The fluorescence resonance energy transfer experiment also showed a decrease in the affinity of glutamine to GlnBP upon mutation of Arg75.24 A mutation at Asp10 or Lys115 disrupted the hydrogen-bond network related to the door keeper and destabilized the ligand interactions in S2LS and S2L, leading the system to S2U or S3U, in which only the ligand interactions at L1,2 were maintained. (In the present 50 ns simulations, these mutants did not open the hinge completely, as RMSDprot < 6 Å, or progressed to S2U directly before reaching S3U.) Finally, the ligand specificity was examined. We performed two additional 50 ns simulations for GlnBP complexed with asparagine (Asn, shorter than glutamine by one methylene group in the side chain), starting from a model structure generated by fitting the carboxylic acid group of Asn to that of the glutamine in the GlnBP bound/closed structure (Figure

Figure 5. (A) Glutamine- and asparagine-bound structures of GlnBP. (B) RMSDprot and the minimum distances for hydrogen bonds in L1,2 (dHB,L1,2), L3,4 (dHB,L3,4), S5,6,7 (dHB,S5,6,7), and #1 and 5 (related to door-keeper formation, dHB,1,5) are plotted as a function of time. Simulation of glutamine-bound GlnBP (gray) and two simulations of asparagine-bound GlnBP (red and blue).

large fluctuations in GlnBP (Figure 5B). The hydrogen-bond network related to the door keeper was almost completely missing during the simulation. This implies that stable ligand binding is tightly coupled to the formation of the door keeper. The length of glutamine is most suitable for binding to both L1,2 and L3,4, but the asparagine side chain is not long enough to form these hydrogen bonds stably. This observation is H

DOI: 10.1021/acs.jpcb.6b11696 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

A. K. was supported by MEXT Grants-in-Aid for Scientific Research (B), 16H04780, and for Challenging Exploratory Research, 16K14714. The computations were performed at the RIKEN Integrated Cluster of Clusters (RICC) and the Research Center for Computational Science, Okazaki, Japan. This research also used the computational resources of the HPCI system, provided by Kyushu University and Nagoya University through the HPCI System Research Project (Project ID: hp140056 and hp150045).

consistent with the weak binding of Asn that was shown in the experimental study.24 The simulations of the complexes with glutamate and aspartate also showed a similar instability of the ligand interaction and protein structure (Figure S7). Therefore, the hydrogen bonds listed in Table 1 ensure not only ligand stability but also ligand specificity.

■

CONCLUSIONS Enhanced sampling of GlnBP complexed with its ligand glutamine was performed at an all-atom resolution by the MSES simulation and revealed the free-energy landscape of the protein−ligand interaction process coupled with the large protein structural change. The free-energy landscape along the reaction coordinates of the protein structure and the protein− ligand interaction described in atomic detail revealed a funnel feature going down from the unbound/open state to the native basin of the bound/closed state. From the classification of the structural ensemble, a set of the substates were determined, with distinctive structural characteristics, and they demonstrated the structural basis of the ligand-binding pathway: first, glutamine associates with Arg75 on the large domain of GlnBP, through strong polar interactions. Subsequently, the ligand interactions with the small domain cause the closing motion of GlnBP, accompanied by desolvation occurring on the protein− ligand interface, which constitutes a rate-limiting step. Finally, the door keeper between the two domains (Asp10−Lys115) and the hydrophobic sandwich between the glutamine ligand and Phe13/Phe50 are formed, to complete the native-specific complex structure. The simulations of the mutants revealed the roles of the ligand-binding sites in stabilizing the ligand interactions. The unstable interactions that appeared in the simulations of GlnBP complexed with asparagine, instead of glutamine, suggested the drastic change in the FES through a small chemical difference in the ligand molecule and the high ligand specificity of GlnBP.

■

■

(1) Shan, Y. B.; Kim, E. T.; Eastwood, M. P.; Dror, R. O.; Seeliger, M. A.; Shaw, D. E. How Does a Drug Molecule Find Its Target Binding Site? J. Am. Chem. Soc. 2011, 133, 9181−9183. (2) Buch, I.; Giorgino, T.; De Fabritiis, G. Complete Reconstruction of an Enzyme-inhibitor Binding Process by Molecular Dynamics Simulations. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 10184−10189. (3) Dror, R. O.; Pan, A. C.; Arlow, D. H.; Borhani, D. W.; Maragakis, P.; Shan, Y. B.; Xu, H. F.; Shaw, D. E. Pathway and Mechanism of Drug Binding to G-protein-coupled Receptors. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 13118−13123. (4) Tiwary, P.; Limongelli, V.; Salvalaglio, M.; Parrinello, M. Kinetics of Protein-ligand unbinding: Predicting Pathways, Rates, and Ratelimiting Steps. Proc. Natl. Acad. Sci. U.S.A. 2015, 112, E386−E391. (5) Limongelli, V.; Marinelli, L.; Cosconati, S.; La Motta, C.; Sartini, S.; Mugnaini, L.; Da Settimo, F.; Novellino, E.; Parrinello, M. Sampling Protein Motion and Solvent Effect During Ligand Binding. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 1467−1472. (6) Held, M.; Metzner, P.; Prinz, J. H.; Noe, F. Mechanisms of Protein-Ligand Association and Its Modulation by Protein Mutations. Biophys. J. 2011, 100, 701−710. (7) Silva, D. A.; Bowman, G. R.; Sosa-Peinado, A.; Huang, X. H. A Role for Both Conformational Selection and Induced Fit in Ligand Binding by the LAO Protein. PloS Comput. Biol. 2011, 7, No. e1002054. (8) Held, M.; Noe, F. Calculating Kinetics and Pathways of Proteinligand Association. Eur. J. Cell Biol. 2012, 91, 357−364. (9) Gu, S.; Silva, D. A.; Meng, L. M.; Yue, A.; Huang, X. H. Quantitatively Characterizing the Ligand Binding Mechanisms of Choline Binding Protein Using Markov State Model Analysis. PloS Comput. Biol. 2014, 10, No. e1003767. (10) Plattner, N.; Noe, F. Protein Conformational Plasticity and Complex Ligand-binding Kinetics Explored by Atomistic Simulations and Markov Models. Nat. Commun. 2015, 6, 7653. (11) Koshland, D. E. Application of a Theory of Enzyme Specificity to Protein Synthesis. Proc. Natl. Acad. Sci. U.S.A. 1958, 44, 98−104. (12) Monod, J.; Wyman, J.; Changeux, J. P. On the Nature of Allosteric Transitions: A Plausible Model. J. Mol. Biol. 1965, 12, 88− 118. (13) Moritsugu, K.; Terada, T.; Kidera, A. Scalable Free Energy Calculation of Proteins via Multiscale Essential Sampling. J. Chem. Phys. 2010, 133, No. 224105. (14) Moritsugu, K.; Terada, T.; Kidera, A. Disorder-to-order Transition of an Intrinsically Disordered Region of Sortase Revealed by Multiscale Enhanced Sampling. J. Am. Chem. Soc. 2012, 134, 7094− 101. (15) Moritsugu, K.; Terada, T.; Kidera, A. Multiscale Enhanced Sampling Driven by Multiple Coarse-grained Models. Chem. Phys. Lett. 2014, 616−617, 20−24. (16) Moritsugu, K.; Terada, T.; Kidera, A. Energy Landscape of AllAtom Protein-Protein Interactions Revealed by Multiscale Enhanced Sampling. PloS Comput. Biol. 2014, 10, No. e1003901. (17) Moritsugu, K.; Terada, T.; Kidera, A. Multiscale Enhanced Sampling for Protein Systems: An Extension via Adiabatic Separation. Chem. Phys. Lett. 2016, 661, 279−283. (18) Fukunishi, H.; Watanabe, O.; Takada, S. On the Hamiltonian Replica Exchange Method for Efficient Sampling of Biomolecular

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b11696. Additional descriptions of the structural representations and simulation data analyses on the time courses and probabilities of the protein−ligand interactions (Table S1 and Figures S1−S7) (PDF)

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +8145-508-7233. Fax: +81-45-508-7367. ORCID

Kei Moritsugu: 0000-0002-5025-6579 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We acknowledge support from the MEXT Creation of Innovation Centers for Advanced Interdisciplinary Research Areas Program in the Project for Developing Innovation Systems and from the Platform for Drug Discovery, Informatics, and Structural Life Science. K. M. was supported by a MEXT Grant-in-Aid for Young Scientists, 15K18520, and I

DOI: 10.1021/acs.jpcb.6b11696 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Systems: Application to Protein Structure Prediction. J. Chem. Phys. 2002, 116, 9058−9067. (19) Sun, Y. J.; Rose, J.; Wang, B. C.; Hsiao, C. D. The Structure of Glutamine-binding Protein Complexed with Glutamine at 1.94 A Resolution: Comparisons with Other Amino Acid Binding Proteins. J. Mol. Biol. 1998, 278, 219−29. (20) Hsiao, C. D.; Sun, Y. J.; Rose, J.; Wang, B. C. The Crystal Structure of Glutamine-Binding Protein from Escherichia Coli. J. Mol. Biol. 1996, 262, 225−242. (21) Oh, B. H.; Pandit, J.; Kang, C. H.; Nikaido, K.; Gokcen, S.; Ames, G. F. L.; Kim, S. H. 3-Dimensional Structures of the Periplasmic Lysine Arginine Ornithine-Binding Protein with and without a Ligand. J. Biol. Chem. 1993, 268, 11348−11355. (22) Oh, B. H.; Kang, C. H.; Debondt, H.; Kim, S. H.; Nikaido, K.; Joshi, A. K.; Ames, G. F. L. The Bacterial Periplasmic HistidineBinding Protein - Structure/Function Analysis of the Ligand-Binding Site and Comparison with Related Proteins. J. Biol. Chem. 1994, 269, 4135−4143. (23) Yao, N.; Trakhanov, S.; Quiocho, F. A. Refined 1.89-Angstrom Structure of the Histidine-Binding Protein Complexed with Histidine and Its Relationship with Many Other Active-Transport Chemosensory Proteins. Biochemistry 1994, 33, 4769−4779. (24) Gruenwald, K.; Holland, J. T.; Stromberg, V.; Ahmad, A.; Watcharakichkorn, D.; Okumoto, S. Visualization of Glutamine Transporter Activities in Living Cells Using Genetically Encoded Glutamine Sensors. Plos One 2012, 7, No. e38591. (25) Chodera, J. D.; Shirts, M. R. Replica Exchange and Expanded Ensemble Simulations as Gibbs Sampling: Simple Improvements for Enhanced Mixing. J. Chem. Phys. 2011, 135, No. 194110. (26) Liu, J. S. Monte Carlo Strategies in Scientific Computing, 2nd ed.; Springer: New York, 2002. (27) Rosso, L.; Tuckerman, M. E. An Adiabatic Molecular Dynamics Method for the Calculation of Free Energy Profiles. Mol. Simul. 2002, 28, 91−112. (28) Maragliano, L.; Vanden-Eijnden, E. A Temperature Accelerated Method for Sampling Free Energy and Determining Reaction Pathways in Rare Events Simulations. Chem. Phys. Lett. 2006, 426, 168−175. (29) Morishita, T.; Itoh, S. G.; Okumura, H.; Mikami, M. Freeenergy Calculation via Mean-force Dynamics Using a Logarithmic Energy Landscape. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2012, 85, No. 066702. (30) 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. (31) Abrams, J. B.; Tuckerman, M. E. Efficient and Direct Generation of Multidimensional Free Energy Surfaces via Adiabatic Dynamics without Coordinate Transformations. J. Phys. Chem. B 2008, 112, 15742−15757. (32) Rosso, L.; Minary, P.; Zhu, Z. W.; Tuckerman, M. E. On the Use of the Adiabatic Molecular Dynamics Technique in the Calculation of Free Energy Profiles. J. Chem. Phys. 2002, 116, 4389−4402. (33) Yamamori, Y.; Kitao, A. MuSTAR MD: Multi-scale Sampling Using Temperature Accelerated and Replica Exchange Molecular Dynamics. J. Chem. Phys. 2013, 139, No. 145105. (34) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-chain Torsion Potentials for the Amber ff99SB Protein Force Field. Proteins 2010, 78, 1950− 1958. (35) Zheng, W.; Brooks, B. R.; Hummer, G. Protein Conformational Transitions Explored by Mixed Elastic Network Models. Proteins 2007, 69, 43−57. (36) Fiser, A.; Sali, A. MODELLER: Generation and Refinement of Homology-based Protein Structure Models. Macromol. Enzymol. 2003, 374, 461−491. (37) The University of Tokyo. http://www.mu2lib.org/ (accessed February 1, 2014).

(38) Jorgensen, W. L.; Chandrasekhar, J; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (39) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-alkanes. J. Comput. Phys. 1977, 23, 327−341. (40) Juffer, A. H.; Botta, E. F. F.; Vankeulen, B. A. M.; Vanderploeg, A.; Berendsen, H. J. C. The Electric Potential of a Macromolecule in a Solvent: a Fundamental Approach. J. Comput. Phys. 1991, 97, 144− 171. (41) 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. (42) Maragliano, L.; Cottone, G.; Ciccotti, G.; Vanden-Eijnden, E. Mapping the Network of Pathways of CO Diffusion in Myoglobin. J. Am. Chem. Soc. 2010, 132, 1010−1017. (43) Moritsugu, K.; Koike, R.; Yamada, K.; Kato, H.; Kidera, A. Motion Tree Delineates Hierarchical Structure of Protein Dynamics Observed in Molecular Dynamics Simulation. Plos One 2015, 10, No. e0131583.

J

DOI: 10.1021/acs.jpcb.6b11696 J. Phys. Chem. B XXXX, XXX, XXX−XXX