Case Study on Temperature-Accelerated Molecular Dynamics

Jun 18, 2014 - This article is part of the International Conference on Theoretical and High Performance Computational Chemistry Symposium special issu...
0 downloads 0 Views 5MB Size
Download PDF
Article pubs.acs.org/JPCA

Case Study on Temperature-Accelerated Molecular Dynamics Simulation of Ligand Dissociation: Inducer Dissociation from the Lac Repressor Protein Yue Hu† and Haiyan Liu*,†,‡,§ †

School of Life Sciences, ‡Hefei National Laboratory for Physical Sciences at the Microscales, and §Hefei Institutes of Physical Science, Chinese Academy of Sciences, University of Science and Technology of China, 96 Jinzhai Road, Hefei, Anhui 230026, China S Supporting Information *

ABSTRACT: We studied ligand dissociation from the inducer-binding domain of the Lac repressor protein using temperature-accelerated molecular dynamics (TAMD) simulations. With TAMD, ligand dissociation could be observed within relatively short simulation time. This allowed many dissociation trajectories to be sampled. Under the adiabatic approximation of TAMD, all but one degree of freedom of the system were sampled from usual canonical ensembles at room temperature. Thus, meaningful statistical analyses could be carried out on the trajectories. A systematic approach was proposed to analyze possible correlations between ligand dissociation and fluctuations of various protein conformational coordinates. These analyses employed relative entropies, allowing both linear and nonlinear correlations to be considered. Applying the simulation and analysis methods to the inducer binding domain of the Lac repressor protein, we found that ligand dissociation from this protein correlated mainly with fluctuations of side-chain conformations of a few residues that surround the binding pocket. In addition, the two binding sites of the dimeric protein were dynamically coupled: occupation of one site by an inducer molecule could significantly reduce or slow down conformational dynamics around the other binding pocket.

1. INTRODUCTION Ligand dissociation often plays important roles in protein function. From the known structure of a protein−ligand complex, computer simulation can be applied to reconstruct possible pathways for ligand dissociation and to analyze necessary conformational changes of the protein along the pathways. However, because of the time-scale limitation of equilibrium molecular dynamics (MD) simulation,1 ligand dissociation usually cannot be directly observed in straightforward simulations. Accelerated sampling methods such as steered MD2 or targeted MD3 have been utilized to study dissociation pathways of biological complexes. For examples, Yang et al. applied steered MD to investigate the disassociation pathway of a type-II kinase inhibitor from its target protein kinase.4 They found that the inhibitor dissociated through a substrate channel rather than an allosteric pocket channel. On the basis of targeted MD and kinetic measurements, Penniston et al. suggested alternative pathways for the association and disassociation between calmodulin and a calmodulin-binding domain.5 In both steered MD and targeted MD, the system is forced to move toward a certain end structure by applying external forces, which are essentially directed to a single target point in a conformational space of many dimensions. The resulting pathway could be distorted by the applied forces. In the meantime, it could be difficult to predict to what extent such distortions might affect simulation results. For examples, © 2014 American Chemical Society

steered MD or targeted MD often lead to unreasonably high energy barriers6,7 (i.e., the calculated work done by the external forces was often much larger than the work expected for an equilibrium or quasi-equilibrium process). Thus, it is useful to explore alternative accelerated simulation approaches to the investigation of dissociation pathways. Unlike steered MD or targeted MD, in temperatureaccelerated MD (TAMD),8−11 the acceleration of sampling is achieved not by artificially changing the potential energy surface of the system but by increasing the temperature. In addition, to avoid unnecessary extensions of the accessible configurational space, not the entire system is simulated at a higher temperature. Instead, one extends the system by adding a set of virtual configurational degrees of freedoms (DOFs), each of which is coupled to a real system DOF that is targeted for acceleration. While the real system is simulated at room temperature, the virtual DOFs are simulated at a much higher temperature, bringing the real DOFs associated with them to cross high free-energy barriers.8 Under the adiabatic decoupling condition,9−11 namely, when the dynamics of the virtual DOFs Special Issue: International Conference on Theoretical and High Performance Computational Chemistry Symposium Received: April 20, 2014 Revised: June 17, 2014 Published: June 18, 2014 9272

dx.doi.org/10.1021/jp503856h | J. Phys. Chem. A 2014, 118, 9272−9279

The Journal of Physical Chemistry A

Article

the relative entropy (RE) or the Kullback−Leibler divergence21−23 between the distributions of each DOF conditioned on different ligand positions. On the basis of such systematic analysis, protein DOFs that significantly correlate with ligand disassociation may be selected, yielding unbiased and detailed pictures about protein dynamics potentially connected to ligand disassociation. We note that the analysis of correlated conformational changes in protein simulations based on relative entropy has been previously proposed by McClendon et al.24,25 From the dissociation trajectories, we observed that the presence or absence of an inducer molecule in the other binding pocket of the protein dimer had a large effect on conformational relaxation in the binding pocket in concern. Namely, there is dynamic coupling between the two binding sites, in accordance with the hypothesis that inducer binding mainly changes conformational plasticity of the protein.

is much slower than that of their associated real DOFs, the virtual DOFs (and the associated real DOFs as well) will sample the high-temperature canonical distribution dictated by the low-temperature potential of mean forces (PMF) of the real system. An important property of TAMD as compared with steered MD or targeted MD is that the underlying dynamics is still driven by the intrinsic potential energy surface undistorted by external forces or restraining potentials. When applied to study dissociation pathways, this property can make TAMD less prone to possible biases caused by preassumed directions of conformational changes. In this work, we applied TAMD to simulate ligand dissociation from the complex between the repressor protein of the lac operon of E. coli (LacR) and its inducer, isopropyl βD-1-thiogalactopyranoside (IPTG).12 In the classical picture of LacR-regulated inducible gene expression, the inducer-free LacR binds to its operator sequence on DNA and represses gene transcription. When bound to inducer molecules such as IPTG, the affinity of LacR to its operator DNA segments is significantly reduced, leaving the repressor protein to drop off from DNA. Then, respective downstream genes can be transcribed.13 According to previous analyses of sequence, structure, and functions,14−17 the entire LacR protein is composed of a N terminal headpiece, followed by an inducerbinding domain comprising a N-terminal and a C-terminal subdomain and then a C-terminal tail. The headpiece is responsible for the binding of operator DNA. The inducerbinding domain is responsible for inducer binding as well as for protein dimerization. The C-terminal tail is responsible for tetramer formation. A range of structural analysis has been carried out on this system, including crystal structures of the inducer-free dimeric protein with the headpiece in complex with operator DNA as well as crystal structures of the dimeric protein without the headpiece and DNA.17 The latter structures include both the inducer-bound and the inducer-free forms. Interestingly, when not in complex with DNA, the inducerbound and the inducer-free structures are almost identical. In the meantime, they all differ significantly from the structure of the DNA-bound form. It was thus proposed that to bind to DNA the protein needs to adopt a conformation that is different from the equilibrium conformation of the DNA-free protein. Inducer binding may not change the equilibrium conformation of the DNA-free protein. Instead, it mainly affects the dynamics or conformational plasticity of the protein,18 making it difficult for the protein to change into the conformation for DNA binding. In the current study, we focused on the process of ligand dissociation from the inducerbinding domain. We simulated the dissociation of one inducer molecule from one monomer site of the dimeric protein with the binding site of the other protein monomer being either empty (noted as the other site empty or OSE) or occupied by another inducer molecule (noted as the other site occupied or OSO). We applied TAMD to sample complete dissociation trajectories of IPTG from LacR. From the trajectories, intermediate ligand positions were clustered using spectral clustering.19 Then, we analyzed possible correlations between protein conformational changes and ligand binding positions along the dissociation pathways. Different conformational DOFs of the protein were examined, including global collective coordinates obtained using principal component analysis20 (PCA), as well as local DOFs describing peptide backbone and side chain conformations. For the latter DOFs, we considered

2. METHODS 2.1. Sampling Dissociation Trajectories by TAMD. The DOF targeted for acceleration by TAMD is the distance d from the center of geometry of the ligand to the center of the ligand binding pocket of the protein. Namely, simulations were performed on an extended system with the following potential energy function 1 Vextended(⇀ r , D) = Vphysical(⇀ r ) + kcoupling[d(⇀ r ) − D]2 2 (1)

r represents coordinates of physical atoms, D in which ⇀ represents a virtual DOF associated with the collective r ), Vphysical and Vextended represent the potential coordinate d(⇀ energy functions of the physical system and of the extended system, respectively, while kcoupling is a force constant that has r ) and D been chosen large enough to keep the values of d(⇀ close to each other. The extended system was simulated with the virtual DOF D sampled at a high temperature of 1600 K by Langevin dynamics, while the physical system was maintained at a temperature of 300 K in normal MD. A dissociation trajectory was sampled by starting a TAMD simulation from a r ) (the initial values of D bound configuration with a small d(⇀ r )) and ending the simulation when d(⇀ r) were set to d(⇀ reached 15 Å. This TAMD approach was applied to sample 100 trajectories of the dissociation of one IPTG molecule from its binding site in a dimeric LacR in explicit solvent. In 50 of the trajectories, the ligand binding site of the other monomer was occupied by another IPTG molecule (the OSO group), while in the other 50 trajectories, the other ligand binding site was empty (the OSE group). Starting configurations for the TAMD simulations have been sampled with a normal equilibrium MD simulation of the dimeric IPTG-LacR complex in solution. (See the simulation details.) The center of the ligand binding pocket has been defined as the center of geometry of the Cα atoms of amino acid residues surrounding the pocket,12 including Leu74, Ala75, Pro76, Ile79, Asn125, Leu148, Asp149, Phe161, Ser193, Arg197, Trp220, Asn246, Asp274, Gln291, Phe293, and Leu296. 2.2. Cluster Ligand Positions along the Dissociation Trajectories. Ligand positions relative to the protein were clustered based on kinetic information contained in the sampled dissociation trajectories, in a spirit similar to the Markov state model26,27 of molecular trajectories. First, the 3-D Cartesian space spanned by the center of the ligand (after 9273

dx.doi.org/10.1021/jp503856h | J. Phys. Chem. A 2014, 118, 9272−9279

The Journal of Physical Chemistry A

Article

different groups based on the clustering of ligand positions, the distribution of a DOF within one group of configurations (i.e., associated with one cluster of ligand position) was compared with the distribution of the same DOF within another group of configurations. We considered the following Jensen−Shannon divergence,22,35,36 which is a symmetrized version of the Kullback−Leibler divergence or relative entropy between two probability distributions

necessary coordinate transformations to superimpose the receptor protein) was divided into small cubic boxes with a side length of 1 Å. (This value has been tested and verified; see the Supporting Information.) Then, configurations contained in the trajectories were grouped into microstates according to which small box the center of the ligand fell into. After that, an intermicrostate transition count matrix S was constructed to record the number of direct transitions from configurations of one group into configurations of another group in the trajectories. The matrix S has been symmetrized by adding together the numbers of transitions in both directions. With the S treated as a similarity matrix and the algorithm of Laplacian spectral clustering19,28−30 applied, the microstates of ligand positions were clustered to obtain a more coarse-grained representation of ligand positions. To be more specific, the transition count matrix S is simplified into a connecting matrix W using the following rules ⎧ 0, if i = j or Sij = 0 ⎪ Wij = ⎨ ⎪ ⎩1, if i ≠ j and Sij > 0

DJS(μ , A , B) =

⎡

pA (μc )

⎢⎣

1/2[pA (μc ) + pB (μc )]

∑ ⎢pA (μc ) ln μC

+ pB (μc ) ln

⎤ ⎥ 1/2[pA (μc ) + pB (μc )] ⎥⎦ pB (μc )

(4)

where A and B represent different groups of conformations with the ligand at different positions, μ represents a DOF describing protein local conformation (e.g., the main chain conformation or side chain conformation of an amino acid residue), and the summation is over all possible states of that DOF, each indicated by μc. The probabilities PA(μc) and PB(μc) were determined from a limited number of sampled configurations and could contain large statistical errors. Thus, the calculated absolute values of DJS(μ, A, B) could not be employed directly to tell if the distribution of the conformational DOF μ changes significantly when the ligand position changes from A to B. We used the following resampling approach to determine the statistical significance of DJS(μ, A, B): the same group of intermediate ligand positions (A or B) is randomly divided into two subgroups (A1 and A2, or B1 and B2) and DJS(μ, A1, A2) or DJS(μ, B1, B2) could be calculated; random division was repeated many times; and the resulting distribution of DJS(μ, A1, A2) and DJS(μ, B1, B2) were compared and the one associated with a larger variation was used as a background relative entropy distribution for the estimation of significance. Finally, the significance or P value associated with DJS(μ, A, B) was calculated as the probability for the relative entropy to be larger than DJS(μ, A, B) in this background distribution. 2.4. Details of the Simulations. Structure of the LacRIPTG complex was obtained from Protein Data Bank (PDB entry 2P9H).17 The CHARMM27 force field37 was used. The CGenFF/ParamChem38 was used to derive force-field parameters for IPTG. (The resulting force-field parameters together with the chemical structure of IPTG are given in Supplementary Tables S1 and S2 and Supplementary Figure S1 in the Supporting Information.) The protonation states of titratable side chains were defined by assuming neutral aqueous environment; namely, the acidic amino acid residues (Asp and Glu) were negatively charged, and the basic amino acid residues (Lys and Arg) were positively charged. Residue His74 was kept positively charged to form the ion pair His74-Asp278 according to experimental results.12 The remaining histidine residues were assigned neutral side chains with proton positions manually determined based on possible hydrogen bonding with surrounding groups. Eight sodium ions were included as counterions to maintain overall neutrality. The TIP3P model was used for water. Water molecules contained in the PDB structure were preserved, and additional water molecules were added to simulate the aqueous environment. Periodic boundary conditions were applied with a rectangular periodic box, with its

(2)

where the i and j are indices of the cubic boxes or microstates. A normalized Laplacian matrix (or the so-called random walk Laplacian matrix19,28) of a graph with the microstates as nodes and Wij as indicators of edges is given by Lrd = I − Δ−1W

1 2

(3)

in which I is the identity matrix and Δ is the diagonal degree matrix of the graph with element Δij = ∑jWijδij, where δij is the Kronecker’s δ symbol. Clustering of the microstates was achieved by diagonalizing Lrd, obtaining its first k eigenvectors with the smallest eigenvalues, followed by treating each microstate as a point in a k-dimensional space with coordinates given by the loads of that microstate in the k eigenvectors. The points (microstates) were then clustered using the k-means algorithm.31 The number of clusters k was chosen naturally based on the relatively large gap between the kth and the (k + 1)th smallest eigenvalues of Lrd. 2.3. Analyzing Correlations between Protein Conformations and Ligand Positions. To describe possible global conformational variations of the protein, we obtained collective coordinates by applying principal component analysis20 (PCA) to the covariant matrix of protein’s Cα atom positions. Configurations from the sampled dissociation trajectories were projected into 2-D principal component spaces to visualize potential differences in protein conformations when the ligand moved to different positions. To examine local conformations of the protein, the main-chain Ramachandran torsional angles as well as the side-chain rotamer states of individual amino acid residues were considered. Main-chain conformations were separated into a number of regions or states on the Ramachandran map (three states for glycine residues, two states for proline residues, and four states for the remaining types of amino acid residues). Discrete rotamer states for side-chain conformations were defined as proposed by Dunbrack et al.32 For side chains containing more than two rotatable torsional angles, only the first two torsional angles were considered. We used the Kullback−Leibler divergence33,34 or relative entropy (RE) to measure if a local conformational DOF significantly correlated with the position of the ligand during the ligand dissociation process. After separating the recorded configurations in the dissociation trajectories into 9274

dx.doi.org/10.1021/jp503856h | J. Phys. Chem. A 2014, 118, 9272−9279

The Journal of Physical Chemistry A

Article

exponential decays (Figure 1B). This behavior suggested that the different TAMD simulations in each group were homogeneous in terms of kinetic rate, and the trajectories were sampled from the same kinetic processes with a single effective rate limiting barrier. Although the time constants of the decays in Figure 1A,B are artificial because of the temperature accelerated dynamics, the decay constants for the OSO and OSE cases are similar, suggesting similar barrier heights. 3.2. Clusters of Ligand Positions. The size of the boxes to discretize ligand positions and the type of Laplacian matrix for spectral clustering was varied, and their effects are discussed in the Supporting Information. On the basis of Supplementary Figures S2 and S3 in the Supporting Information, a box side length of 1 Å and the random walk Laplacian were chosen for final analysis. On the basis of an obvious gap between the eighth and the ninth eigenvalues of the Laplacian matrix for clustering (see Supplementary Figure S3C in the Supporting Information), ligand positions along the disassociation trajectories were grouped into eight clusters. (See Supplementary Figure S3D in the Supporting Information.) The ligand positions belonging to different clusters are indicated together with the crystal structure of the protein in Figure 2.

initial side lengths determined by setting the minimum distance from protein atoms to box boundaries to 12 Å. Smooth Particle-Mesh Ewald (PME) method on GPU39 with grid sizes of 96 × 96 × 96 was used to compute the electrostatic interactions. The PME calculations were performed every two integration steps. A cutoff of 9 Å and a cutoff of 7.5 Å of a switch function were applied to handle the van der Waals interactions. All covalent bonds formed by a hydrogen atom with a heavy atom were constrained at their ideal lengths. Masses of hydrogen atoms were multiplied by a factor of four so that a relatively long integration time step of 4 fs could be used. Masses of the heavy atoms to which the hydrogen atoms were attached were reduced accordingly so that the total mass of the system was unchanged. After setup, the system was subjected to 500 steps of energy minimization, followed by 100 ps of NVT MD simulation in which the system was gradually heated to 300 K and then another 1 ns NPT simulation for equilibration. After that, a 200 ns NPT simulation was performed to sample initial structures for subsequent TAMD sampling. The 50 initial structures for the other-site-occupied or OSO simulations were extracted evenly from the sampled time trajectory. To start the other-siteempty or OSE group of TAMD simulations, one IPTG was removed from the extracted structures. The TAMD simulations lasted a few nanoseconds to 100 ns before the ligand could escape 15 Å away from the native binding sites. Configurations were recorded every 4 ps for subsequent analysis. MD simulations were carried out with the program ACEMD,40 enabling efficient uses of GPU hardware. Selfimplemented codes were used to perform TAMD every time step using the plugin interface provided by ACEMD (source codes of the plugin could be obtained from the authors on request). The PYMOL41,42 or the VMD43 programs were used for molecular graphics. VMD, Carma,44 or Wordom45,46 were used for other analyses.

Figure 2. Clustered ligand positions (small crosses) relative to the crystal structure of the inducer-binding domain of LacR in the TAMD trajectories. Ligand positions in different clusters are in different colors: cluster 1 in black, cluster 2 in orange, cluster 3 in blue, cluster 4 in red, cluster 5 in green, cluster 6 in yellow, cluster 7 in gray, and cluster 8 in purple. (A,B) Views of the same structure in different orientations, with panel B rotated from panel A by 90° around a horizontal axis. Small red spheres indicate the locations of IPTG in the crystal structure.

3. RESULTS AND DISCUSSION 3.1. Dissociation Trajectories. With TAMD, the ligand escaped from the receptor in a few to 100 ns, while in the 200 ns normal MD simulations the ligand always stayed inside the binding pocket. Thus, high temperature for one DOF was able to accelerate the dissociation process sufficiently. Figure 1A shows the numbers of remaining trajectories (i.e., trajectories in which the ligand remained within 15 Å from the center of the binding pockets) in the OSO and OSE groups, respectively, as a function of simulation time. The curves fit well to simple

Among the clusters of ligand positions, cluster 7 corresponds to the native binding position. The other clusters deviated from this position in two directions: one went toward the front side of the crystal structure shown in Figure 2, namely, from cluster 7 to cluster 3 and then clusters 1 and 8; the other went toward the backside of the shown structure, that is, from cluster 7 to cluster 4 and then to clusters 2 and 5. Dissociation in the frontward direction can be regarded as dominating, with cluster 3 being much larger than cluster 4. Indeed, in 97 of the total 100 sampled trajectories, the ligand dissociated along the frontward pathway. Four movies showing example dissociation processes along the two pathways in the two groups of simulations are provided in the Supporting Information. In Figure 3, representative structures with ligand positions belonging to cluster 7 (the native binding position), cluster 3 (the frontward intermediate position), and cluster 4 (the backward intermediate position) were shown. For comparisons, structures from the OSO and OSE groups of trajectories were

Figure 1. (A) Numbers of trajectories with ligand remained within 15 Å from the center of the inducer binding pocket of LacR as a function of the length of the TAMD simulations. (B) Same as panel A, but logarithms of the numbers and their fitting straight lines are shown. OSO or OSE stand for TAMD simulations with the other inducer binding site occupied (OSO) or empty (OSE). 9275

dx.doi.org/10.1021/jp503856h | J. Phys. Chem. A 2014, 118, 9272−9279

The Journal of Physical Chemistry A

Article

Figure 3. Representative structures showing the dissociating inducer at different positions relative to the protein environment. Panels A, D, and G show structures in cluster 7 (the native binding position). Panels B, E, and H show structures in cluster 3 (the frontward intermediate position). Panels C, F, and I show structures in cluster 4 (the backward intermediate position). In each of panels A−C, two structures, one from the OSO group (magenta) and another from the OSE group (cyan) of simulations, were superimposed. Panels D−I show the same structures as panels A−C, respectively, with enlarged views, surrounding residues in contact with the ligand within 4 Å shown as sticks. Panels D−F show structures from the OSO group. Panels G−I show structures from the OSE group. Atoms are colored by element type: O in red, C in green, N in blue, and S in orange. Hydrogen atoms are not shown for clarity.

tions of JS divergences by repetitively and randomly splitting the sample set associated with cluster 3 into two sets and calculating the JS divergences between the two sets. (Similar splitting of the sample set associated with cluster 7 was also performed, and the resulting JS divergences had a much narrower distribution.) In Figure 4A,B, the JS divergences for the backbone distributions between cluster 7 and cluster 3 are compared with the estimated background JS divergences for the OSO and OSE trajectories, respectively. From these Figures, we could not find any particular residues whose backbone conformations were strongly correlated with ligand position. In Figure 4C,D, JS divergences for side-chain conformations are compared. Interestingly, the JS divergences associated with only a few residues exhibited significantly larger values than corresponding background averages. The differences can be confirmed by P values (Figure 4E,F). The residues included Asp149 and Trp220 for the OSO group and Asp149 for the OSE group. The locations of the above residues in the 3-D structure of the protein were inspected. Asp149 can be considered as on the inner wall of the inducer binding pocket (Figures 3). Its side chains directly form hydrogen bonds with the inducer molecule. Trp220 is located at the entrance of the binding pocket (Figure 5). Its rotation seems to be necessary for the inducer to move between the bound position cluster 7 and the intermediate position cluster 3. It is thus no coincidence that

superimposed in Figure 3A−C. The residues in direct contacts with the ligand at different positions are shown in remaining subfigures (Figure 3D−I). The OSO and OSE groups of trajectories were quite similar in terms of their intermediate ligand positions and the group of protein residues interacting with the ligand at these positions. 3.3. Possible Correlations between Ligand Position and Protein Conformation. The cluster 7 (the native binding position) and cluster 3 (the intermediate position for frontward disassociation) are large enough (Supplementary Table S3 in the Supporting Information) to allow meaningful analysis of possible correlations between protein conformation and ligand position. First, possible correlations involving global conformational degrees of freedom (DOFs) of the protein were examined by principle component analysis (PCA). The protein global motions show no significant correlations with ligand position by this method (Supplementary Figure S4 in the Supporting Information). To look for possible correlations between ligand position and protein local conformation, we considered relative entropies or JS divergences between respective distributions of local backbone or side-chain conformational states. Again, the two clusters of ligand position compared were cluster 7 (the native binding position) and cluster 3 (the intermediate position for frontward dissociation). To consider significance of the computed JS divergences, we obtained background distribu9276

dx.doi.org/10.1021/jp503856h | J. Phys. Chem. A 2014, 118, 9272−9279

The Journal of Physical Chemistry A

Article

affinity and exhibited decreased allosteric amplitude.48 In ref 49, Gatti-Lafranconi et al. investigated the mutant Trp220Phe and found the mutant repressor had significantly reduced leakiness, namely, more stable binding with DNA in the absence of an inducer.49 Thus, the mutations of Asp149 or Trp220 had significant impacts on the phenotype of LacR. We would like to note that such impacts cannot be connected with the TAMD analysis directly, although the TAMD results do suggest these amino acid residues to be among the more important ones for inducer-binding related functions. Another interesting note to which we will come back in discussions is that the correlations between ligand position and local side-chain conformations are statistically much more significant in the OSO group than in the OSE group of trajectories, the differences being qualitative in the cases of Trp220.

4. DISCUSSION The systematic analyses of the TAMD dissociation trajectories suggested that inducer dissociation from the inducer−protein complex was not associated with large global conformational changes of the dimeric protein. This was the case for both the OSO and OSE groups of trajectories. This observation is consistent with previous structural studies showing that the inducer-bound and inducer-free proteins have similar structures when not bound to DNA.17 In terms of protein local conformational DOFs correlated with ligand disassociation, the OSO and OSE groups of trajectories exhibited subtle but potentially important differences. In the OSO groups of simulations, the side-chain conformation distribution of Trp220 had significant correlations with ligand positions. A more detailed examination of the distributions suggested that in the OSO group of trajectories the side-chain conformation of Trp220 was less dynamic and more restricted when the ligand was within cluster 3 than it was within cluster 7 (Supplementary Figure S5A in the Supporting Information). The correlation diminished in the OSE groups of simulations mainly because such restriction did not seem to exist there (Supplementary Figure S5B in the Supporting Information). This suggests differences in conformational dynamics around the ligand binding site between the OSO and OSE groups of trajectories. To further look into this issue independent of the somewhat artificial clustering of ligand positions, we inspected how the distributions of the distance from the center of the ligand to the center of the binding pocket evolved with time in the TAMD simulations. Two sets of time windows were considered to investigate the time evolutions. The first set was based on the absolute simulation time: each recorded configuration frame was assigned to one of seven time windows based on the absolute time at which the frame was sampled. Each of the first six time windows covered a 10 ns section, and the seventh time window covered the remaining time, which was >60 ns. By this definition of time windows, coordinate frames in shorter-lived trajectories contributed to the statistics of only the earlier but not the later time windows. The second set of time windows was defined based on relative time: the sampling time for each recorded frame was scaled by the lifetime of the respective trajectory. The resulting relative time, with a value range from 0 to 1, was divided evenly into seven windows, and each recorded configuration frame was assigned to one of the windows accordingly. In Figure 6A,B, the distance distributions in different absolute time windows were compared between the

Figure 4. (A) For the OSO group of trajectories, the JS divergences (red lines) between the distributions of backbone conformational state with the ligand position in cluster 7 (the native binding position) and with the ligand position in cluster 3 (the intermediate position for frontward dissociation). The averages and standard variations of the background JS divergences are shown in blue and in green, respectively. (B) Same as panel A but for the OSE trajectories. (C) Same as panel A but for distributions of side-chain conformations. (D) Same as panel C but for the OSE trajectories. (E,F) P values (cyan lines) computed from the distributions in panels C and D, respectively, and the corresponding JS divergences (red lines).

Figure 5. Trp220 is located at the entrance of the inducer binding pocket. (A) Two monomers, one in space-filling model and the other in ribbons, of the LacR dimer are shown. Trp220 in the first monomer is in red color. (B) Trp220 (sticks) in the crystal structure (green) and in a representative structure with ligand position in cluster 7 (the native binding position) (yellow).

the side-chain conformations of these residues were singled out to have significant correlations with ligand dissociation. There were in fact experimental reports of the large effects of site-directed mutations of the above residues on inducer binding.47−49 Suckow et al. reported phenotypes of 4000 mutants of LacR.47 The amino acid residues at positions 149 and 220 both fell into a group of residues whose mutations are likely to lead to the so-called Is phenotype; namely, the mutants are no longer responsive to inducer. In ref 48, Xu et al. created mutants at either or both positions 125 and 149. The mutants were found to bind to IPTG with five to eight times lower 9277

dx.doi.org/10.1021/jp503856h | J. Phys. Chem. A 2014, 118, 9272−9279

The Journal of Physical Chemistry A

Article

In summary, through the case study on IPTG dissociation from LacR, we have shown that TAMD can be an efficient method to sample dissociation trajectories without applying external forces that might distort the potential energy surface. From the TAMD trajectories, correlations between ligand positions and protein conformational dynamics can be reliably analyzed. Finally, we note that there are still some important limitations of the current study, some of them awaiting future work. For example, from the current TAMD trajectories, we cannot extrapolate a free-energy profile because the simulations did not sample the state in which the bound and unbound forms were in equilibrium. In addition, the current TAMD strategy may not be suitable to sample the association process because unlike dissociation association is an entropy-decreasing process and may not be easily accelerated by raising the temperature.

■

ASSOCIATED CONTENT

S Supporting Information *

Structure, topology, and parameters of IPTG. Effect of grid size on the clustering of ligand positions. Effects of using different Laplacian matrices on clustering results. Cluster sizes. Possible correlations between ligand position and protein global conformation. Distributions of side chain torsional angles χ1 and χ2 of Trp220. Example dissociation trajectories, with filenames indicating trajectory group and direction of ligand movement. This material is available free of charge via the Internet at http://pubs.acs.org.

Figure 6. (A) Distributions of the distance between the centers of the ligand and the binding pocket within different absolute time windows in the OSO trajectories. (B) Same as panel A but for the OSE trajectories. (C) Same as panel A but for different relative time windows. (D) Same as panel C but for the OSE trajectories. (E) Distributions of the radius of gyration of residues surrounding the inducer binding pocket in different relative time windows in the OSO trajectories. (F) Same as panel E but for the OSE trajectories.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: 86-551-3607450. Fax: 86-5513607451. Notes

The authors declare no competing financial interest.

OSO and OSE groups of trajectories. Figure 6C,D compared similar distributions in different relative time windows. Described with either set of time windows, the time evolution of the OSO distributions showed shifts of peak positions accompanied by peak broadening, while the time evolution of the OSE distributions showed no shift of peak positions but only peak broadening. This indicated that in the OSO case, the free-energy surface underlying the ligand dissociation process was less smooth than in the OSE case. Figure 6E,F demonstrates the same point with the time evolution of another geometric coordinate, namely, the radius of gyration of the amino acid residues surrounding the binding pocket. This coordinate has been examined because it could reflect the “effective size” of the binding pocket. In both the OSO and the OSE trajectories, the distributions did not change much with time. However, the OSO and OSE distributions were different in that the OSO distributions showed a clear two-peak structure, while the minor peak in the OSE distributions was almost flat, suggesting a smoother free-energy surface for protein local conformational fluctuations in the OSE case. The difference between the OSO and OSE trajectories should be attributed to the presence or absence of an inducer molecule at the other binding site. The previous results suggest that there is dynamic coupling between the two binding sites of the dimeric protein: protein conformational dynamics around one binding site from which ligand dissociation was examined may be significantly reduced or slowed down by the occupation of the other binding site by an inducer.

■

ACKNOWLEDGMENTS This work has been supported by Chinese Natural Science Foundation (grant number 21173203).

■

REFERENCES

(1) van Gunsteren, W. F.; Bakowies, D.; Baron, R.; Chandrasekhar, I.; Christen, M.; Daura, X.; Gee, P.; Geerke, D. P.; Glattli, A.; Hunenberger, P. H.; et al. Biomolecular Modeling: Goals, Problems, Perspectives. Angew. Chem., Int. Ed. 2006, 45 (25), 4064−4092. (2) Isralewitz, B.; Baudry, J.; Gullingsrud, J.; Kosztin, D.; Schulten, K. Steered Molecular Dynamics Investigations of Protein Function. J. Mol. Graphics Modell. 2001, 19 (1), 13−25. (3) Schlitter, J.; Engels, M.; Kruger, P.; Jacoby, E.; Wollmer, A. Targeted Molecular-Dynamics Simulation of Conformational Change Application to the T↔R Transition in Insulin. Mol. Simul. 1993, 10 (2−6), 291−308. (4) Yang, L. J.; Zou, J.; Xie, H. Z.; Li, L. L.; Wei, Y. Q.; Yang, S. Y. Steered Molecular Dynamics Simulations Reveal the Likelier Dissociation Pathway of Imatinib from Its Targeting Kinases C-Kit and Abl. PLoS One 2009, 4 (12), e8470. (5) Penniston, J. T.; Caride, A. J.; Strehler, E. E. Alternative Pathways for Association and Dissociation of the Calmodulin-Binding Domain of Plasma Membrane Ca(2+)-Atpase Isoform 4b (Pmca4b). J. Biol. Chem. 2012, 287 (35), 29664−29671. (6) Yang, K.; Liu, X.; Wang, X.; Jiang, H. A Steered Molecular Dynamics Method with Adaptive Direction Adjustments. Biochem. Biophys. Res. Commun. 2009, 379 (2), 494−498.

9278

dx.doi.org/10.1021/jp503856h | J. Phys. Chem. A 2014, 118, 9272−9279

The Journal of Physical Chemistry A

Article

(30) Ng, A. Y.; Jordan, M. I.; Weiss, Y. On Spectral Clustering: Analysis and an Algorithm. In Advances in Neural Information Processing Systems 14, Vols. 1 and 2; MIT Press: Cambridge, MA, 2002; Vol. 14, pp 849−856. (31) Hartigan, J. A.; Wong, M. A. Algorithm as 136: A K-Means Clustering Algorithm. Appl. Stat. 1979, 100−108. (32) Dunbrack, R. L., Jr. Rotamer Libraries in the 21st Century. Curr. Opin. Struct. Biol. 2002, 12 (4), 431−440. (33) Kullback, S. Information Theory and Statistics; Courier Dover Publications: New York, 2012. (34) Kullback, S.; Leibler, R. A. On Information and Sufficiency. Ann. Math. Stat. 1951, 79−86. (35) Lin, J. H. Divergence Measures Based on the Shannon Entropy. IEEE Trans. Inf. Theory 1991, 37 (1), 145−151. (36) Menéndez, M.; Pardo, J.; Pardo, L.; Pardo, M. The JensenShannon Divergence. J. Franklin Inst. 1997, 334 (2), 307−318. (37) 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 (18), 3586−3616. (38) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. Charmm General Force Field: A Force Field for Drug-Like Molecules Compatible with the Charmm All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31 (4), 671−690. (39) Harvey, M. J.; De Fabritiis, G. An Implementation of the Smooth Particle Mesh Ewald Method on Gpu Hardware. J. Chem. Theory Comput. 2009, 5 (9), 2371−2377. (40) Harvey, M. J.; Giupponi, G.; De Fabritiis, G. Acemd: Accelerating Biomolecular Dynamics in the Microsecond Time Scale. J. Chem. Theory Comput. 2009, 5 (6), 1632−1639. (41) DeLano, W. L. Pymol Molecular Viewer: Updates and Refinements. Abstr. Pap. Am. Chem. Soc. 2009, 238. (42) DeLano, W. L.; Lam, J. W. Pymol: A Communications Tool for Computational Models. Abstr. Pap. Am. Chem. Soc. 2005, 230, U1371−U1372. (43) Humphrey, W.; Dalke, A.; Schulten, K. Vmd: Visual Molecular Dynamics. J. Mol. Graphics Modell. 1996, 14 (1), 33−38. (44) Glykos, N. M. Software News and Updates - Carma: A Molecular Dynamics Analysis Program. J. Comput. Chem. 2006, 27 (14), 1765−1768. (45) Seeber, M.; Cecchini, M.; Rao, F.; Settanni, G.; Caflisch, A. Wordom: A Program for Efficient Analysis of Molecular Dynamics Simulations. Bioinformatics 2007, 23 (19), 2625−2627. (46) Seeber, M.; Felline, A.; Raimondi, F.; Muff, S.; Friedman, R.; Rao, F.; Caflisch, A.; Fanelli, F. Wordom: A User-Friendly Program for the Analysis of Molecular Structures, Trajectories, and Free Energy Surfaces. J. Comput. Chem. 2011, 32 (6), 1183−1194. (47) Suckow, J.; Markiewicz, P.; Kleina, L. G.; Miller, J.; KistersWoike, B.; MullerHill, B. Genetic Studies of the Lac Repressor 0.15. 4000 Single Amino Acid Substitutions and Analysis of the Resulting Phenotypes on the Basis of the Protein Structure. J. Mol. Biol. 1996, 261 (4), 509−523. (48) Xu, J.; Liu, S.; Chen, M. Z.; Ma, J. P.; Matthews, K. S. Altering Residues N125 and D149 Impacts Sugar Effector Binding and Allosteric Parameters in Escherichia Coli Lactose Repressor. Biochemistry 2011, 50 (42), 9002−9013. (49) Gatti-Lafranconi, P.; Dijkman, W. P.; Devenish, S. R. A.; Hollfelder, F. A Single Mutation in the Core Domain of the Lac Repressor Reduces Leakiness. Microb. Cell Fact. 2013, 12, 67−77.

(7) Noe, F.; Fischer, S. Transition Networks for Modeling the Kinetics of Conformational Change in Macromolecules. Curr. Opin. Struct. Biol. 2008, 18 (2), 154−162. (8) 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 (1−3), 168−175. (9) 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 (11), 4389−4402. (10) VandeVondele, J.; Rothlisberger, U. Canonical Adiabatic Free Energy Sampling (Cafes): A Novel Method for the Exploration of Free Energy Surfaces. J. Phys. Chem. B 2002, 106 (1), 203−208. (11) Hu, Y.; Hong, W.; Shi, Y. Y.; Liu, H. Y. TemperatureAccelerated Sampling and Amplified Collective Motion with Adiabatic Reweighting to Obtain Canonical Distributions and Ensemble Averages. J. Chem. Theory Comput. 2012, 8 (10), 3777−3792. (12) Lewis, M. The Lac Repressor. C. R. Biol. 2005, 328 (6), 521− 548. (13) Jacob, F.; Monod, J. Genetic Regulatory Mechanisms in the Synthesis of Proteins. J. Mol. Biol. 1961, 3 (3), 318−356. (14) Bell, C. E.; Lewis, M. A Closer View of the Conformation of the Lac Repressor Bound to Operator. Nat. Struct. Biol. 2000, 7 (3), 209− 214. (15) Bell, C. E.; Lewis, M. The Lac Repressor: A Second Generation of Structural and Functional Studies. Curr. Opin. Struct. Biol. 2001, 11 (1), 19−25. (16) Bell, C. E.; Lewis, M. Crystallographic Analysis of Lac Repressor Bound to Natural Operator O1. J. Mol. Biol. 2001, 312 (5), 921−926. (17) Daber, R.; Stayrook, S.; Rosenberg, A.; Lewis, M. Structural Analysis of Lac Repressor Bound to Allosteric Effectors. J. Mol. Biol. 2007, 370 (4), 609−619. (18) Hawkins, R. J.; McLeish, T. C. Coarse-Grained Model of Entropic Allostery. Phys. Rev. Lett. 2004, 93 (9), 098104. (19) von Luxburg, U. A Tutorial on Spectral Clustering. Stat Comput 2007, 17 (4), 395−416. (20) Jolliffe, I. Principal Component Analysis; Wiley Online Library: Encyclopedia of Statistics in Behavioral Science: New York, 2005; Vol. 3, pp 1580−1584. (21) Kullback, S. The Kullback-Leibler Distance. In American Statistician; American Statistical Association: Washington, DC, 1987; Vol. 41, pp 340−340. (22) Lindorff-Larsen, K.; Ferkinghoff-Borg, J. Similarity Measures for Protein Ensembles. PLoS One 2009, 4, 1. (23) Qian, H. Relative Entropy: Free Energy Associated with Equilibrium Fluctuations and Nonequilibrium Deviations. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2001, 63 (4), 042103. (24) McClendon, C. L.; Friedland, G.; Mobley, D. L.; Amirkhani, H.; Jacobson, M. P. Quantifying Correlations between Allosteric Sites in Thermodynamic Ensembles. J. Chem. Theory Comput. 2009, 5 (9), 2486−2502. (25) McClendon, C. L.; Hua, L.; Barreiro, G.; Jacobson, M. P. Comparing Conformational Ensembles Using the Kullback-Leibler Divergence Expansion. J. Chem. Theory Comput. 2012, 8 (6), 2115− 2126. (26) Pande, V. S.; Beauchamp, K.; Bowman, G. R. Everything You Wanted to Know About Markov State Models but Were Afraid to Ask. Methods 2010, 52 (1), 99−105. (27) Bowman, G. R.; Beauchamp, K. A.; Boxer, G.; Pande, V. S. Progress and Challenges in the Automated Construction of Markov State Models for Full Protein Systems. J. Chem. Phys. 2009, 131 (12), 124101. (28) Meila, M.; Shi, J. B. Learning Segmentation by Random Walks. Advances in Neural Information Processing Systems 13 2001, 13, 873− 879. (29) Shi, J. B.; Malik, J. Normalized Cuts and Image Segmentation. Ieee Transactions on Pattern Analysis and Machine Intelligence 2000, 22 (8), 888−905. 9279

dx.doi.org/10.1021/jp503856h | J. Phys. Chem. A 2014, 118, 9272−9279