Protein–Ligand Dissociation Simulated by Parallel Cascade

This is an open access article published under a Creative Commons Attribution (CC-BY) License, which permits unrestricted use, distribution and reproduction in any medium, provided the author and source are cited.

Article Cite This: J. Chem. Theory Comput. 2018, 14, 404−417

pubs.acs.org/JCTC

Protein−Ligand Dissociation Simulated by Parallel Cascade Selection Molecular Dynamics Duy Phuoc Tran,† Kazuhiro Takemura,‡ Kazuo Kuwata,§ and Akio Kitao*,∥ †

Downloaded via 141.101.132.239 on July 6, 2018 at 13:15:47 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

Department of Computational Biology and Medical Sciences, Graduate School of Frontier Sciences, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa-shi, Chiba 277-8562, Japan ‡ Institute of Molecular and Cellular Biosciences, The University of Tokyo, 1-1-1 Yayoi, Bunkyo-ku, Tokyo 113-0032, Japan § Center for Emerging Infectious Diseases, Gifu University, 1-1 Yanagido, Gifu-shi, Gifu 501-1194, Japan ∥ School of Life Science and Technology, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8550, Japan S Supporting Information *

ABSTRACT: We investigated the dissociation process of tri-N-acetyl-D-glucosamine from hen egg white lysozyme using parallel cascade selection molecular dynamics (PaCS-MD), which comprises cycles of multiple unbiased MD simulations using a selection of MD snapshots as the initial structures for the next cycle. Dissociation was significantly accelerated by PaCS-MD, in which the probability of rare event occurrence toward dissociation was enhanced by the selection and rerandomization of the initial velocities. Although this complex was stable during 1 μs of conventional MD, PaCS-MD easily induced dissociation within 100−101 ns. We found that velocity rerandomization enhances the dissociation of triNAG from the bound state, whereas diffusion plays a more important role in the unbound state. We calculated the dissociation free energy by analyzing all PaCS-MD trajectories using the Markov state model (MSM), compared the results to those obtained by combinations of PaCS-MD and umbrella sampling (US), steered MD (SMD) and US, and SMD and the Jarzynski equality, and experimentally determined binding free energy. PaCS-MD/MSM yielded results most comparable to the experimentally determined binding free energy, independent of simulation parameter variations, and also gave the lowest standard errors. molecular dynamics (SMD) with US12−14 utilize MD for calculating free energy differences. Sufficient sampling at a reasonable computational cost can be achieved by using enhanced sampling techniques, for example, replica exchange molecular dynamics (REMD),15 metadynamics,8 simulated annealing methods,16 and parallel cascade selection molecular dynamics (PaCS-MD).17, These and other methods help the system avoid being trapped in local free energy minima or metastable states, freeing them to sample new conformations or states. Combinations of SMD with US (SMD/US),14 MD with restraining potentials,18−21 replica exchange umbrella sampling (REUS),22 and targeted MD (TMD)23 with US (TMD/US)24 can simulate ligand dissociation from proteins and calculate the binding free energy at a reasonable computational cost. In SMD/US, a bias force is applied to generate ligand dissociation pathways in SMD, and US samples local conformational spaces divided into multiple windows along the pathway. The binding free energy is obtained as the potential of mean force between the bound and unbound states. For the complex formed between lysozyme and its antibody (HyHEL-10), the protein structures were distorted when the bias force was applied to the center of mass (COM) of the proteins, whereas the obtained

1. INTRODUCTION Calculation of free energy differences between distinct molecular states has a long history and remains an active research theme in computational chemistry, physics, and biophysics. Computational methods for determining free energy can be categorized into four major groups: thermodynamic integration (TI), sampling-based methods, nonequilibrium dynamics, and adaptive biasing techniques.1 Kirkwood2 first introduced TI, which is widely used. TI makes use of the adiabatic evolution of statistical average along reaction coordinates to retrieve free energy differences. Zwanzig proposed alchemical free energy perturbation, which decomposes a free energy change into multiple intermediate steps.3 Later, the free energy perturbation approach was extended to methods based on reaction coordinates, including umbrella sampling (US).4,5 Of the nonequilibrium dynamics methods, the Jarzynski equality6 is the most widely used to determine the relation between free energy change and nonequilibrium trajectories. Adaptive biasing dynamics methods, e.g., metadynamics7,8 and Wang−Landau methods,9 monitor the reaction coordinates and penalize the visited region by using adaptive forces.1 Molecular dynamics (MD) simulation is a valuable method for calculating the free energy. For example, free energy perturbation,10 thermodynamics integration,11 US with the weighted histogram analysis method (WHAM),4 and steered © 2017 American Chemical Society

Received: May 15, 2017 Published: November 28, 2017 404

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Article

Journal of Chemical Theory and Computation

and MSM allows a more cost-effective and accurate evaluation of the binding free energy.

potential of mean force (PMF) was improved if large distortion was carefully avoided in the multistep TMD.24 These results indicate that careful treatment is required to reduce artifacts due to force biases, thus motivating us to conduct dissociation simulation without applying force biases. In this work, we employed PaCS-MD17,25 to generate ligand dissociation pathways without applying force biases. PaCS-MD incorporates cycles of selection of initial structures for multiple independent MD simulations and conformational sampling with MD. The initial structures for MD simulation in each cycle are selected based on the ranking of an appropriate quantity that characterizes conformational transitions so that the expected conformation gradually approaches the target. The selection of snapshots enhances conformational changes by increasing the probability of rare event occurrences to drive the conformation toward the target. PaCS-MD17 was originally designed to sample conformational space toward specified target structures. The sampled conformational transition pathways were then further analyzed by additional US along the pathways. Nontargeted PaCS-MD (nt-PaCS-MD) enhances sampling without specifying the target by selecting the most deviated structures from the average structure based on GramSchmidt orthogonalization.25 The outlier flooding method (OFLOOD)26 can also be considered an extension of PaCSMD and uses outlier detection based on clustering. Temperature-aided cascade molecular dynamics enhances conformational sampling with independent MD simulations at different temperatures.27 Peng and Zhang developed the PaCS-Fit method to conduct selection based on experimental data and successfully applied it to fit small-angle X-ray scattering and electron microscopy data.28 In this paper, we demonstrate that PaCS-MD can be used to simulate protein−ligand dissociation within tens of nanoseconds by employing a longer intermolecular distance as the target quantity for the selection of the initial structures without applying force bias. PaCS-MD-generated dissociation pathways are compared to those from SMD. We utilized all the PaCSMD trajectories to calculate free energy change along the dissociation pathways in combination with the Markov state model (MSM). For comparison, free energy calculations were also performed using combinations of PaCS-MD and US, SMD and US, and SMD and the Jarzynski equality. We examined the dissociation of tri-N-acetyl-D-glucosamine (triNAG) from hen egg white lysozyme (LYZ). LYZ has long been the subject of many studies due to its antibacterial property, as recently described.29 triNAG binds to a cleft between two domains: a domain consisting of α helices (α domain) and a β-rich domain (β domain). Experimental and computational studies showed that the cleft has six N-acetyl-Dglucosamine (NAG) binding pockets, designated A-F, and that binding to A-B-C is the main binding mode.30−33 Recently, Zhong and Pastel used a polarizable force field, together with molecular mechanics and generalized Born and surface area (MM-GBSA), to investigate the A-B-C and B-C-D binding modes of triNAG to LYZ. However, neither of their models reproduced the binding free energy of the wild-type LYZtriNAG complex,34 which suggests the need to include the contribution of the solvation free energy to the binding free energy. We show that the main interactions between LYZ and triNAG in the bound state agree with those in the crystal structure and that the obtained binding free energy of the LYZtriNAG complex is in agreement with experimental and other computational results. Moreover, the combination of PaCS-MD

2. MATERIALS AND METHODS 2.1. Parallel Cascade Selection Molecular Dynamics. PaCS-MD generates conformational transition pathways using cycles of multiple and/or parallel short conventional MD simulations.17 The simulation starts with multiple independent MD simulations without bias. The simulated systems in these MDs (hereinafter called “replicas”) are identical except for the initial velocities and conformations. The number of MD simulations conducted is nrep, and the MD simulation time for each MD is tcyc. The initial structures for the next cycle are the top nrep structures selected from the ranking of the generated snapshots based on a certain quantity. Typically, this is a position along the reaction coordinates. The new MD cycle is restarted after rerandomization of the initial velocities with the Maxwell−Boltzmann distribution. The cycle of MD and selection is repeated until the generated snapshots closely approach the target value. We used the inter-COM distance between LYZ and triNAG to rank the snapshots for the selection. Since the LYZ-triNAG complex is in the global free energy minimum conformation, movement toward dissociation rarely occurs. Dissociating motion is enhanced by selecting snapshots with longer interCOM distances. The use of PaCS-MD should allow the generation of more natural dissociation pathways for triNAG compared to SMD, because the former method does not apply any bias to the system. We considered the initial snapshots as well as the MD-generated snapshots in each cycle as candidates for the next cycle. Therefore, if no increase in the inter-COM distance was observed, then the same initial snapshots were employed for the next cycle but with different initial velocities. One trial of PaCS-MD thus generates a set of unbiased trajectories that highly overlap in conformational space along the triNAG dissociation pathway, since new replica MDs are started from snapshots of the previous cycle. Each trajectory gives information on a sampled local conformational space. We employed MSM to combine the information on overlapping local conformational spaces and to build consistent information on overall conformational space along the dissociation pathway. We also generated joint trajectories connecting the initial bound state and the final unbound state along the dissociation pathways, which are termed “reactive trajectories” in the original PaCS-MD paper (see Figure 1b of ref 17). A reactive trajectory is generated by joining the fragments of MD trajectories that start from the initial snapshots and end at the selected snapshots. Therefore, the number of reactive trajectories obtained in one trial of PaCS-MD is equal to nrep. The reactive trajectories were used as representative dissociation pathways and are also employed as initial structures for US. The reactive trajectories can also be ordered based on the ranking of the last cycle. 2.2. Steered Molecular Dynamics. In addition to PaCSMD, SMD35 was also conducted to dissociate triNAG from LYZ. SMD is an in silico method for stretching a molecule from its complex, inspired by single-molecule force spectroscopy such as atomic force microscopy (AFM) or optical tweezer experiments.36 During SMD, LYZ was weakly restrained around its COM, and a dummy atom, which moved with constant velocity along the z axis, is connected to the COM of triNAG via a virtual spring. This type of SMD is called the constant velocity pulling method. The dummy atom can be 405

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Article

Journal of Chemical Theory and Computation

Figure 1. Visualization of the simulation box in the initial state and the key amino acid residues of LYZ that interact with triNAG after equilibration. (a) Overall arrangement of the LYZ-triNAG complex and solvent in the simulation box and (b) a close-up view. The residues shown by yellow Licorice models are disulfide-bonded cysteine residues, and the molecule represented as a multicolored Licorice model is triNAG. (c) A view along the z-axis to show the electrostatic potential on the LYZ surface (blue: positive charges, red: negative charges). triNAG is shown as a Licorice model. (d) LigPlot+ diagram to show the interactions between LYZ and triNAG. Hydrophobic contacts are represented as spline curves outlining residue labels, and hydrogen bonds are shown as dotted lines. (e) Positions of the LYZ residues involved in hydrogen bonds with triNAG in the binding pockets. Blue and red residue labels show the residues situated in the α and β domains, respectively. Panels (a-c,e) and (d) were created by VMD77 and LigPlot+,78 respectively.

considered as the tip in AFM experiments, and its initial position is the same as COM of triNAG. The spring is 1 represented by a harmonic potential, 2 k[ξ(t ) − λ(t )]2 , where k is a force constant, ξ(t) is the inter-COM distance between LYZ and triNAG at SMD time t, λ(t) = λ0 + vt is the distance between the COM of LYZ and the dummy atom (λ0 is the initial distance), and v is the pulling velocity corresponding to the constant velocity of the dummy atom. SMD trajectories were used to calculate the free energy profile of the dissociation in combination with the Jarzynski equality and were also employed as the initial structures for US. 2.3. Markov State Model. The MSM introduces a discrete-state stochastic kinetic model of the observed process and is a powerful tool to obtain insights for linking experimental and simulation data.37,38 There are three steps to build an MSM from MD trajectories: building microstates, building the transition matrix, and validating the built MSM. To build the microstates, high dimensional data is projected onto lower dimensional space by clustering, principle component analysis (PCA),39 and time-lagged independent component analysis (TICA).40−42 In this study, we applied K-means clustering43 for the inter-COM distance for building the microstates, and then the transition matrix between the microstates

T = {Tij} = {P[x(t + τ ) ∈ Sj|x(t ) ∈ Si]}

(1)

was estimated. Each matrix element Tij is calculated between a pair of microstates (i and j) at a predetermined lag time τ, using the maximum likelihood estimation procedure described in ref 44. In eq 1, x(t) and x(t + τ) represent the coordinates at time t and (t + τ). We examined the independence of the implied time scale from τ and chose a suitable value of τ for building the MSM. The equilibrium probabilities of the microstates can be calculated as the stationary eigenvector p = {pi} of T, and the equilibrium free energy of microstate i can be obtained as −(ln pi)/β.45 We used this method to calculate the free energy profile as a function of the inter-COM distance. To build each MSM, we employed all the full MD trajectories generated by each trial of PaCS-MD, regardless of whether or not the snapshots of the trajectories are selected for the next cycle. It should be noted that the selected and nonselected trajectories together provide significant information to estimate transition probabilities between microstates. We calculated the inter-COM distance between LYZ and triNAG in each snapshot in all trajectories of a trial to build a single data set for MSM. This data set is used as input for MSM construction as mentioned in the above procedure. From the built MSMs of all trials, the obtained equilibrium probabilities as a function of the inter-COM distance are averaged and shown as a distribution of free energy along the inter-COM 406

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Article

Journal of Chemical Theory and Computation Table 1. List of Simulations and Their Conditions simulation

nrep/tcyc (ns) in PaCS-MD or v in SMD (nm/ns)

no. of trials

10/0.1 100/0.1 10/1.0 1.25 0.25 0.05

10 10 10 24 12 12

10,0.1

PaCS-MD PaCS-MD100,0.1 PaCS-MD10,1 SMDfast SMDmed SMDslow

MD time to reach 4 nm (ns)a 3.5 1.5 24.6 2.0 9.8 45.9

± ± ± ± ± ±

1.0 0.2 10.1 0.1 0.4 9.3

MSMb

USb

Jarzynskib

Y (5.42) Y (2.00) Y (2.67)

Y (7.67) Y (12.12) Y (6.28) Y (6.96)

Y (0.12) Y (0.28) Y (0.96)

a MD simulation time per replica averaged over trials. The values after “±” show the standard deviation. b‘Y’ indicates that dissociation free energy profile was calculated with a combination of the specified simulation method and free energy calculation method. The values in parentheses indicate total simulation time including all the MD trials and sampling simulations in μs.

and NaCl to ensure an ion concentration of 0.15 M and charge neutrality. We used the AMBER99SB force field51 for LYZ and the GLYCAM06 force field52 for triNAG. All simulations were performed by GROMACS 5.0.5.53 The equilibration simulation was conducted using the following steps: 1) Energy minimization was performed with the steepest descent method followed by the conjugate gradient method, imposing heavy atom positional restraints with a force constant of 1000 kJ/mol nm2. 2) The system was heated from 0 K to 300 K in 500 ps using the NVT ensemble and then equilibrated for an additional 500 ps at 300 K. The time step of all the MD simulations was 1.0 fs. 3) MD was conducted using the NPT ensemble for the next 100 ps at 300 K with a relaxation time of 0.1 ps for heat bath coupling and at 1.0 atm with a 2.0 ps relaxation time for isotropic pressure coupling. 4) The force constant of the positional restraints was deducted by 100 kJ/mol nm2 every 100 ps until it vanished. 5) The system was equilibrated for 3.0 ns using the NPT ensemble with weak positional restraints imposed on the sulfur atoms of the cysteine residues (shown in yellow in Figure 1(a)) so that the overall translation and orientation of LYZ was maintained. The LINCS method54 to constrain the bond lengths and the leapfrog integration method55 were used in steps 2−4, while the velocity Verlet method56 without bond constraints was employed in step 5. An isothermal condition was realized by velocity rescaling57 in steps 2 and 3 and using a Nosé−Hoover method58,59 in steps 4 and 5. An isobaric condition was achieved using a Berendsen barostat60 in step 3, a ParinelloRahman barostat61 in step 4, and a MTTK barostat62 in step 5. 2.7. Simulations Performed. The results of step 5 provided the start structures for PaCS-MD and SMD, with the initial structures for PaCS-MD selected from the middle 1 ns trajectory of step 5. The simulations conducted are listed in Table 1. The PaCS-MD cycles were repeated until the shortest inter-COM distances in the top nrep rank exceeded 7 nm. We conducted 10 trials of PaCS-MD with nrep = 10 or 100 and ts = 0.1 ns (PaCS-MD10,0.1 and PaCS-MD100,0.1, respectively) as well as 10 trials with nrep = 10 and ts = 1 ns (PaCS-MD10,1). The initial structures for SMD were selected from the last 2 ns trajectory of step 5. The snapshots were sampled every 0.5 ps, including the initial structures, and therefore 201 and 2001 snapshots from each 0.1 and 1.0 ns trajectory were considered for selection, respectively. SMD was performed with k = 500 kJ/mol nm2 and three different pulling velocities: v = 1.25 (SMDfast), 0.25 (SMDmed), and 0.05 nm/ns (SMDslow). The pulling force was applied along the z-axis. We performed 24 runs for SMDfast and 12 runs each for SMDmed and SMDslow. For US along the reaction coordinates, we selected 50 snapshots from the SMD or PaCS-MD trajectories. Starting

distance reaction coordinate. We used MSMBuilder 3 for MSM46 construction and analysis. 2.4. Umbrella Sampling Using the Weighted Histogram Analysis Method. The free energy profile of LYZtriNAG dissociation was also calculated by US, starting from the snapshots generated by PaCS-MD and SMD (namely PaCS-MD/US and SMD/US, respectively). US calculates free energy from a probability distribution in equilibrium. Restraint MD simulations with the umbrella potential V are conducted around different points along a reaction coordinate. In this work, V was applied to the inter-COM distance between LYZ and triNAG along the dissociation pathway. The free energy profile ΔG(d) can be calculated by WHAM.4 We selected snapshots from the reactive trajectories of PaCS-MD and the trajectories of SMD. The program distributed by the Grossfield Lab47 was used for WHAM. 2.5. Jarzynski Equality. In this work, the free energy profile along the inter-COM distance between LYZ and triNAG was also calculated with the Jarzynski equality from SMD trajectories (SMD/Jarzynski hereafter). The Jarzynski equality6 relates the free energy difference ΔG with the total work performed on the system, W as e−β ΔG = ⟨e−βW ⟩

(2)

where ⟨···⟩ indicates the statistical average. External work can be easily obtained via the following equation:48 W0 → t = −kv

∫0

t

dt ′[ξ(t ′) − λ(t ′)]

(3)

Equations 2 and 3 contribute to the fundamental exact relations linking the irreversible work on nonequilibrium processes to the equilibrium statistics of a given system.49 The second order cumulant expansion was used to calculate PMF from the Jarzynski equality.50 2.6. Simulation Setup. To generate the initial model, we set up a simulation box containing wild-type LYZ and triNAG (PDB ID: 1HEW). The initial box was 7.9094 × 7.66642 × 16.26561 nm3 along the x- y-, and z-axes, respectively, and allows large dissociation movement of triNAG along the z-axis (Figure 1a). Initially, the inter-COM distance between LYZ and triNAG was directed parallel to the z-axis. To avoid significant overall translation and reorientation of LYZ during triNAG dissociation, weak positional restraints were applied to the sulfur atoms of the cysteine residues involved in the four disulfide bonds of LYZ during the final stage of equilibration (step 5; see next paragraph) and in the production runs. As shown in Figure 1(a), this system size was chosen so that the distances between the outermost atoms of the complex and the box edges were at least 1.5, 1.5, and 5.5 nm along the x-, y-, and z-axes, respectively. The box was solvated with TIP3P water 407

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Article

Journal of Chemical Theory and Computation

highest inter-COM distance between LYZ and triNAG during the PaCS-MD simulation is shown in Figure 2. LYZ and

from each snapshot, we conducted 10 ns restrained MD with a force constant of 500 kJ/mol nm2. We ensured the stability of the LYZ-triNAG model throughout long simulation times by performing 1 μs conventional MD simulation at 300 K and 1.0 atm. We also conducted a conventional MD simulation of triNAG in 0.15 M NaCl solution with the TIP3P water model to calculate the selfdiffusion constant of free triNAG and used the same simulation box size to avoid box size effects on diffusion.63−65 In this paper, values after ± show standard deviations unless otherwise specified.

3. RESULTS 3.1. Interactions between LYZ and triNAG in the Bound State and the Stability of the Complex. triNAG consists of three NAG units (NAG I-III) and binds to the A-BC pockets of LYZ. The LYZ residues interacting with triNAG were examined after the relaxation MD (step 5) (Figure 1(c) and 1(d)). The binding cleft of LYZ includes positively charged residues (ARG73 and ARG112), negatively charged residues (ASP52 and ASP101), polar residues (GLN57, ASN59, and ASN103), and hydrophobic residues (LEU56, ILE58, TRP62, TRP63, LEU75, ILE98, ALA107, TRP108, and VAL109). Negative charges (shown in red in Figure 1(c)) are mostly located in the bottom of the binding cleft, while positive charges (blue in Figure 1(c)) are widely distributed on the surface of LYZ, except for the bottom of the cleft surface. Therefore, the negatively charged residues in the cleft are surrounded by the positively charged and hydrophobic residues. triNAG formed long-lasting hydrogen bonds with ASN59, TRP62, TRP63, ASP101, ASN103, and ALA107 of LYZ during the 1 μs conventional MD run, identical to those found in the crystal structure (PDB ID: 1HEW)31 and by other computational studies.32−34 Of the above-listed residues, ASP52 is known to play several key roles in catalysis.66 Furthermore, the isomerization of ASP101 is involved in the aging process, resulting in decreased affinity of LYZ for the triNAG ligand.67 Interestingly, a LYZ mutant in which cavities in the hydrophobic core are filled by swapping amino acids between residues 12 and 56 (MET12 and LEU56 in the wild-type) exhibited increased stability.68 Using multicanonical MD, Kamiya et al. showed that ASP101 and TRP62 are important amino acids for the interaction of triNAG and LYZ during the association process.69 We examined the stability of the LYZ-triNAG complex during a 1 μs conventional MD simulation. We calculated the inter-COM distance, d, as well as the root-mean-square deviation (RMSD) of LYZ and triNAG as a function of time after performing least-squares fitting of the LYZ backbone atoms to the crystal structure (Figure S1). The inter-COM distance between LYZ and triNAG was stable, with a value around 1.41 ± 0.03 nm, although the RMSDs of both LYZ and triNAG slightly shifted after 0.7 μs. The number of hydrogen bonds between LYZ and triNAG was also stable (4.2 ± 1.1) during the simulation. The largest root-mean-square fluctuation (RMSF) of the LYZ atoms from the average was 0.49 nm (for ARG68, located in a loop), and all residues with an RMSF larger than 0.2 nm were located outside the binding cleft. The results of the 1 μs MD simulation clearly showed the stability of the LYZ-triNAG complex during this time scale. 3.2. LYZ-triNAG Dissociation by PaCS-MD. Although the LYZ-triNAG complex was stable during 1 μs, PaCS-MD can dissociate the complex very easily. The time evolution of the

Figure 2. Evolution of the inter-COM distance between lysozyme and triNAG, d, in the top reactive trajectories during each PaCS-MD trial as a function of the number of cycles for (a) PaCS-MD10,0.1, (b) PaCSMD100,0.1, and (c) PaCS-MD10,1. The meanings of the shaded regions are described in the main text.

triNAG completely dissociated (d > 4 nm) at 34.8 ± 10.3 (3.48 ns), 14.5 ± 1.7 (1.45 ns), and 24.6 ± 10.1 cycles (24.6 ns) on average, and the simulations were stopped at 41.6 ± 10.4, 20.2 ± 2.5, and 27.6 ± 10.5 cycles for PaCS-MD10,0.1, PaCSMD100,0.1, and PaCS-MD10,1, respectively, when d reached 7 nm (Figure 2 and Table 1). Compared to PaCS-MD10,0.1, the number of cycles required for complete dissociation was reduced to 48.6 and 66.3% in PaCS-MD100,0.1 and PaCSMD10,1, respectively. PaCS-MD100,0.1 and PaCS-MD10,1 required the same computational resource per cycle, but the sampling efficiency of PaCS-MD100,0.1 was higher than that of PaCSMD10,1, because it required fewer cycles to achieve complete dissociation. The difference in cycle reduction rates indicates that increasing the number of replicas (nrep) is a more efficient way to reduce the number of cycles required to achieve dissociation than is simulation length (tcyc). The standard deviation of the number of cycles required for dissociation is also the smallest for PaCS-MD100,0.1, which minimized the variation between simulation lengths required for multiple trials, as noted in Figure 2. Since the probability of observing rare events is proportional to nrep, the increase in nrep reduces the number of cycles required for dissociation. We classified the dissociation process into three states. The first is the bound state (the regions below the shaded regions in Figure 2), in which the inter-COM distance increases slowly and almost linearly. The second is the partially bound state (the shaded regions in Figure 2: 1.79−3.65, 1.73−3.39, and 1.82− 3.18 nm for PaCS-MD10,0.1, PaCS-MD100,0.1, and PaCS-MD10,1, respectively), where a nonlinear rapid increase of d was observed but triNAG was in partial contact with LYZ. The third is the unbound state (the regions above the shaded regions), in which d increased almost linearly and rapidly. The average total number of cycles required for complete dissociation strongly depends on the number of cycles spent in the bound state, which were 24.1 ± 11.2 (2.41 ns), 7.8 ± 1.8 (0.78 ns), and 17.7 ± 9.5 (17.7 ns) for PaCS-MD10,0.1, PaCS-MD100,0.1, and PaCS408

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Article

Journal of Chemical Theory and Computation MD10,1, respectively. Cycles with no significant increase of d (namely, “trapped” cycles) occurred 8.8 ± 8.4 and 6.6 ± 6.34 times on average in PaCS-MD 10,0.1 and PaCS-MD10,1 , respectively, whereas no trapped cycle was observed in PaCSMD100,0.1. This indicates that PaCS-MD, with a larger nrep, is a more efficient way of avoiding trapped cycles. Traps mostly occurred in the bound and partially bound states. The average number of continuous trapped cycles was 4.4 ± 2.3 and 3.8 ± 2.9 for PaCS-MD10,0.1 and PaCS-MD10,1, respectively. In the unbound state, the average speeds of triNAG movement were 0.64 ± 0.08, 0.77 ± 0.16, and 1.56 ± 0.27 nm/cycle for PaCSMD10,0.1, PaCS-MD100,0.1, and PaCS-MD10,1, which correspond to 6.4, 7.7, and 1.6 nm/ns, respectively. The speed of triNAG movement in the PaCS-MD10,1 simulation was equivalent to the pulling velocity of SMDfast (1.25 nm/ns), while those of PaCSMD10,0.1 and PaCS-MD100,0.1 were 5 times faster than the pulling velocity of SMDfast. 3.3. Effects of Velocity Rerandomization and Selection on triNAG Dissociation during PaCS-MD. We examined diffusive properties during PaCS-MD by inspecting the self-diffusion constant D by the Einstein relation: D = lim t →∞

⟨Δr 2(t )⟩ 6t

COM distance and calculated the self-diffusion constants in each cycle (Dbound, Dpartial, and Dunbound), as depicted in Figure 3.

(4)

D was determined by least-square fitting of the time evolution of mean square displacement (MSD) Δr2(t) to a straight line. We first calculated the effective diffusion in “reactive trajectories”, which were used for the initial structures for US. As mentioned in Section 2.1, the reactive trajectory is a joint multiple trajectory connecting the trajectories from the initial bound state to the final unbound state,17 and the initial structures for PaCS-MD/US were selected from the top reactive trajectories. To analyze the characteristics of the bound, partially bound, and unbound states, each reactive trajectory was split into ranges of the three states based on the inter-COM distance, as shown in Figure 2. The obtained Δr2(t) is shown in Figure S2. Self-diffusion constants in the unbound −5 state, Dreact. unbound, were 6.6 ± 1.9, 7.7 ± 1.2, and 1.9 ± 0.8 × 10 2 10,0.1 100,0.1 10,1 cm /s for PaCS-MD , PaCS-MD , and PaCS-MD , respectively, which indicates that shorter tcyc (= 0.1 ns) accelerated effective diffusion more than 3-fold compared to tcyc = 1 ns. For comparison with the diffusion constants in the unbound state (Dreact. unbound), the diffusion constant of the LYZ-free state was also calculated. The values of Dreact. unbound obtained from PaCS-MD simulations were significantly larger than the obtained free diffusion constant, 1.1 ± 0.5 × 10−5 cm2/s, indicating that PaCS-MD enhanced the effective diffusion constants of the unbound state. Velocity rerandomization causes perturbation of the reactive trajectories at each concatenating point in the trajectory, because they are continuous in conformational space but discontinuous in phase space. We considered another property of the reactive trajectories, namely, the average time length of trajectory fragments in the reactive trajectories, Δtfrag. If Δtfrag is too short, the system might not relax sufficiently after velocity rerandomization. However, the Δtfrag values for PaCS-MD10,0.1 and PaCS-MD10,1 were 79.2 ± 7.7 and 842.4 ± 122.8 ps, respectively, which are significantly longer than the safe limit exchange time interval 4 ps in temperature REMD.70 Next, we analyzed the dynamics of triNAG in individual MD trajectories. We classified the PaCS-MD trajectories into the bound, partially bound, and unbound states based on the inter-

Figure 3. Distributions of the triNAG self-diffusion constants, Dbound (red), Dpartial (blue), and Dunbound (green) in (a) PaCS-MD10,0.1, (b) PaCS-MD100,0.1, (c, d) PaCS-MD10,1 shown as probability densities. The densities in (c) and (d) were calculated from entire 1 ns and the first 0.1 ns trajectories, respectively. Insets show each COM trajectory of triNAG around LYZ (white cartoon model) in different colors, depending on the values of the diffusion constant: blue (D ≤ 0.5; all units 10−5 cm2/s), green (0.5 < D ≤ 1.0), yellow (1.0 < D ≤ 1.5), orange (1.5 < D ≤ 2.0), and red (2.0 < D). The values after ± show standard deviations.

We also show these constants calculated for the first 0.1 ns of the PaCS-MD10,1 trajectories (Figure 3(d)) for comparison with the results of PaCS-MD10,0.1. The distributions of the diffusion constants are very similar to those from simulated random walk trajectories generated by varying the concentration of random point obstacles.71 The values of Dunbound (1.2 ± 0.2, 1.0 ± 0.6, and 1.2 ± 0.7 × 10−5 cm2/s for PaCS-MD10,0.1, PaCS-MD100,0.1, and PaCS-MD10,1) are in good agreement with the free diffusion constant (1.1 ± 0.5 × 10−5 cm2/s). If velocity rerandomization has a significant effect on diffusion, the diffusion coefficient should depend on the simulation length (1.0 or 0.1 ns), because the effect weakens as the MD simulation becomes longer. Therefore, the effect of velocity rerandomization is weak in the unbound state, and no clear artifact was observed in the diffusion constants. In the bound state, Dbound is almost the same for the first 0.1 ns (Figure 3 (a,b,d)) but is smaller for 1.0 ns (Figure 3(c)), which indicates that the effect of velocity rerandomization on diffusion depends on the length of the trajectory in the bound state. Interestingly, Dbound and Dpartial (red and blue curves in Figure 3) are mainly populated below 0.5 × 10−5 cm2/s and spatially form a low mobility region around the binding pockets (the trajectories shown by blue in the insets of Figure 3). As triNAG dissociates further, higher mobility regions were found but were not homogeneously distributed, as clearly seen by color variations. Since the range of Dunbound values is broader than that of the 409

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Article

Journal of Chemical Theory and Computation

near the end of the MD run tended to be selected more frequently, because the movement of triNAG in the unbound state is largely determined by diffusion, and larger deviations should occur near the end of the MD simulation in a diffusiondominant environment. Overall, the results indicate that dynamics in individual short simulations behaved as expected in unbiased MD simulations, and thus we judged that MSM can be appropriately applied to PaCS-MD trajectories. 3.4. LYZ-triNAG Dissociation by SMD. In SMD, ligand dissociation was induced by the steering force. Figure 5 shows

other states, trajectories with different mobilities are mixed in the unbound state. Dunbound of PaCS-MD10,0.1, PaCS-MD100,0.1, and PaCS-MD10,1 is significantly larger than Dbound, by 5.8, 5.3, and 20.5 times, respectively. Now let us consider the effects of MD length, velocity rerandomization, and selection on PaCS-MD sampling. The Dbound and Dunbound values obtained from PaCS-MD10,0.1 (Figure 3(a)) and the first 0.1 ns of PaCS-MD10,1 (Figure 3(d)) agreed very well, which confirms that the behavior of the system in the first 0.1 ns of PaCS-MD10,1 and PaCS-MD10,0.1 simulation was the same. The smaller value of Dbound in PaCS-MD10,1 obtained from full 1.0 ns trajectories (Figure 3(c)) compared to that obtained from the first 0.1 ns indicates that a longer MD simulation time did not accelerate diffusion in the bound state. These results suggest that velocity rerandomization enhanced sampling in the bound state. To study the effect of velocity rerandomization on selection, we calculated the probability of selection as a function of simulation time of the selected snapshots (Figure 4). In the bound state, snapshots near the

Figure 5. (a-f) Time evolution of the inter-COM distance, d, and (g-l) force between LYZ and triNAG as a function of SMD time for (a,d,g,j) SMDfast, (b,e,h,k) SMDmed, and (c,f,i,l) SMDslow. (a-c) and (g-i) show the cases in which single force peaks as a function of time were observed (also see Table 2). (d-f) and (j-l) show the cases in which double force peaks as a function of time were observed. In these cases, small plateau regions are seen in (d-f).

Figure 4. Probability of selection as a fraction of time of the selected snapshots versus the total length of each MD (1.0 or 0.1 ns) in each state for (a) PaCS-MD10,0.1, (b) PaCS-MD100,0.1, and (c) PaCS-MD10,1 in the bound (red), partially bound (blue), and unbound (green) states.

the time evolution of the inter-COM distance, d, and the force between LYZ and triNAG as a function of SMD time. Unlike PaCS-MD, the time evolution of d showed a steep jump between two linear regions (Figure 5(a-f)). The initial linear regions range from 1.3 to 1.7 nm in all three cases. The jump started when the steering force became maximum. The average maximum forces during the dissociation process were 336.1 ± 44.5, 303.9 ± 43.6, and 267.2 ± 33.5 kJ/mol·nm in the SMDfast, SMDmed, and SMDslow simulations, respectively. These values are in the same range as the AFM disruption forces previously reported.72 The lower the pulling rate, the weaker the maximum force required to dissociate triNAG. After the steep jump, d linearly increased, and the force converged toward zero at around 2.2−2.6 nm. We found two different patterns in the steering force as a function of d: single peak and double peaks. Table 2 shows the average values of the inter-COM distance between LYZ-triNAG at the first and second force peaks, the

beginning of each MD tended to be selected, indicating that velocity rerandomization enhances movements leading toward dissociation in the PaCS-MD scheme, but this effect quickly decays. Similar tendencies were observed in the partially bound states during PaCS-MD10,0.1 and PaCS-MD10,1 but not during PaCS-MD100,0.1, likely due to the aforementioned trapped cycle observed in the bound and partially bound states of the structures during PaCS-MD10,0.1 and PaCS-MD10,1. If a significant increase in d was not observed, snapshots near the beginning of the MD run were selected, which raised the probability of snapshot selection from this time region. Since the trap did not occur in PaCS-MD100,0.1, the probability did not peak near the beginning. In the unbound state, snapshots 410

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Article

Journal of Chemical Theory and Computation Table 2. Characteristic Inter-COM Distances and Forces in SMDd 1st peakc sim. SMD

fast

SMDmed SMDslow

a

b

type

no.

single double single double single double

13 11(1) 5 7(1) 8 4(0)

d (nm) 1.5 1.5 1.5 1.5 1.5 1.4

± ± ± ± ± ±

0.1 0.0 0.1 0.1 0.1 0.0

2nd peakc

force (kJ/mol nm) 351.8 321.9 309.4 292.7 269.4 247.9

± ± ± ± ± ±

38.7 40.0 46.6 41.8 37.2 11.5

convergencec

d (nm)

force (kJ/mol nm)

2.2 ± 0.3

182.1 ± 35.1

2.4 ± 0.5

159.4 ± 24.9

2.0 ± 0.3

149.7 ± 5.9

d (nm) 2.6 2.6 2.5 2.4 2.2 2.3

± ± ± ± ± ±

0.2 0.2 0.1 0.1 0.2 0.3

a

single: single peak was seen, double: two peaks were observed. bThe number in parentheses is the number of cases in which the heights of the two peaks are almost equal. cAverage values are shown, and the values after ± show standard deviations. dAverage values of the inter-COM distance between LYZ/triNAG at the first and second force peaks as a function of SMD time, values of force at the peaks, and the distances where the force converged.

values of the force at the peaks, and the distances where the force converged. Note that values in parentheses in Table 2 indicate the number of cases in which the heights of the two peaks are the same (so-called same-height double peaks). We found that the heights of the first and second peaks decreased as the SMD velocity decreased. The force peak for the singlepeak cases was generally larger than that of the double-peak cases. The heights of the same-height double peaks were lower than the single peak by 100 kJ/mol·nm. After reaching the first peak, triNAG quickly dissociated from LYZ; however, triNAG remained trapped in the double-peak cases, which correspond to small plateau regions (red lines in Figure 5). The standard deviations of the positions of the first peaks were always small (≤0.1 nm), showing that the first stage of the dissociation processes in SMD started from the same position. 3.5. Dissociation Pathways in PaCS-MD and SMD. We analyzed the spatial distribution of the dissociation pathways to obtain better insight into the relation between the dissociation pathway and free energy. Figure 6(a) shows the COM positions of triNAG along 10 representative reactive trajectories, each of which is the top ranked reactive trajectory in each PaCSMD10,0.1 trial. The inset of Figure 6(a) depicts the triNAG COM positions in all PaCS-MD trajectories generated in one trial. Since PaCS-MD trajectories are started from selected snapshots of the previous cycle, the generated trajectories mutually overlap in conformational space. Therefore, a set of COM trajectories generated in one trial shows a dissociation pathway which connects the bound and completely unbound states, as shown in the inset of Figure 6(a). We quantified the sample conformational space along each dissociation pathway by calculating an effective diameter σeach(d) of a cross section of trajectories as a function of the inter-COM distance d. σeach(d) is defined as σ each(d) =

4S each(d) π

Figure 6. (a) Dissociation pathways of triNAG represented by the COM positions of triNAG (small spheres) in the first reactive trajectories of 10 PaCS-MD10,0.1 trials from LYZ (white cartoon model). Inset shows all the trajectories generated in a representative trial of PaCS-MD10,0.1 (red in main panel). (b) Effective diameter σ of the sampled area per trial of PaCS-MD as a function of the inter-COM distance d. (c) Effective diameter σ over all trials. The inset shows a close-up. PaCS-MD10,0.1 (red), PaCS-MD100,0.1 (blue), PaCS-MD10,1 (green), SMDfast (magenta), SMDmed (orange), and SMDslow (black). Error bars show standard deviations.

(5)

where Seach(d) is a cross section of the sampled triNAG COM positions at d. The average of σeach(d) over trials is shown in Figure 6(b). σ each (d) is larger in the bound state (d < 2 nm) but is almost flat after complete dissociation (d > 2.5 nm). This reflects that more cycles were spent in the bound state than in the unbound state, resulting in larger σ each (d) values in the bound state. The plateau value of σ̅ in the unbound state is consistent with triNAG diffusing essentially freely in the unbound state. In the unbound state, the average values of σ each (d) were 0.96 ± 0.33, 1.50 ± 0.74, and 2.07 ± 1.31 nm for

PaCS-MD10,0.1, PaCS-MD100,0.1, and PaCS-MD10,1, respectively. This shows that a longer simulation time for each replica provides a larger sampling diameter compared to increasing the number of replicas. In the partially bound state, the sampling 411

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Article

Journal of Chemical Theory and Computation diameter is comparable between PaCS-MD100,0.1 and PaCSMD10,1. We also examined the variation of the dissociation pathways generated by distinct PaCS-MD and SMD trials (Figure S3). This figure clearly shows that significantly different dissociation pathways are generated in each type of simulation. To quantify this variation, we also calculated the σ(d) for the PaCS-MD reactive trajectories of all trials and all SMD trajectories (Figure 6(a)), denoted as σall(d), and the results are shown in Figure 6(c). In the bound state, the average values of σall(d) were 1.03 ± 0.08 and 1.06 ± 0.15 (nm) in the PaCS-MD10,0.1 and PaCSMD100,0.1 simulations, respectively, and are comparable to the values 0.91 ± 0.17 and 1.03 ± 0.41 (nm) obtained by SMDfast and SMDmed, respectively. However, σall(d) obtained by SMDslow and PaCS-MD10,1 for the bound state were significantly larger, 1.76 ± 0.29 and 2.61 ± 0.17 (nm), respectively. In the unbound state, σall(d) obtained by SMDmed and SMDslow steeply increased as d increased. We note that diffusion governs the movement along the x and y directions, because a pulling force was applied only along the z direction. Therefore, the ratio of SMDslow simulation time spent at d = 4 nm versus SMDmed is 4.7, consistent with the ratio of σall(d = 4), 4.4 (Figure 6(c) and Table 1). 3.6. Dissociation Free Energy. The free energy profile (potential of mean force, PMF) of triNAG dissociation from LYZ as a function of the inter-COM distance was calculated by combinations of PaCS-MD and MSM (PaCS-MD/MSM), PaCS-MD and US (PaCS-MD/US), SMD and US (SMD/US), and SMD and the Jarzynski equation (SMD/Jarzynski) (Figure 7). Since the free energy profiles were obtained as the average over distinct dissociation pathways (shown in Figure S3), they should be clearly distinguished from the minimum free energy path. In this paper, we are mainly interested in the free energy difference between two flat regions in the profiles, i.e., the bound state and completely unbound state. The free energy profiles were all stable in the inter-COM distance range 4.0−4.5 nm, and we defined the dissociation free energy ΔGd as the energy difference between the bound state and the average value over d = 4.0−4.5 nm. We expect that the calculated dissociation free energies are equal to the negative value of the binding free energy ΔGb as ΔGb = −ΔGd. In PaCS-MD/MSM, a MSM was constructed using PaCSMD trajectories generated by PaCS-MD100,0.1 or PaCS-MD10,1 simulations (Table 1). Note that the MSM was built using all trajectories of each trial, and the average PMF was obtained over all trials, as shown in Figure 7(a). Each trial of PaCSMD10,0.1 lacked adequate statistics to build the MSM properly. After carefully checking the evolution of the number of microstates and the implied time scale as a function of lag time τ, we determined 50 microstates for both cases and selected 45 and 305 ps as the best τ values for PaCS-MD100,0.1 and PaCSMD10,1, respectively. The number of microstates we tried ranged from 10 to 100. We judged 50 is the best, because distances between neighboring microstates approximately are 0.1 nm, which is suitable to calculate smooth PMF. Also, we obtained enough statistics to calculate transition probabilities. The lag times were chosen to fulfill the requirement of statistics: the minimum of the half of trajectory length and the maximum of lag time before disconnecting the state involved in the slowest process.45 We tried different lag times with 5 ps separation and monitored the results. We confirmed that the PMF curve converges as the lag time approaches the selected value. These lag time values were much shorter than values

Figure 7. Dissociation free energy profiles calculated by combinations of (a) PaCS-MD and MSM (blue: PaCS-MD100,0.1/MSM, green: PaCS-MD10,0.1/MSM), (b) PaCS-MD and US (red: PaCS-MD10,0.1/ US, green: PaCS-MD10,1/US), (c) SMD and US (magenta: SMDfast/ US, orange: SMDmed/US, black: SMDslow/US), and (d) SMD and the Jarzynski equality (magenta: SMDfast/Jarzynski, orange: SMDmed/ Jarzynski, black: SMDslow/Jarzynski) as functions of the inter-COM distance d. Error bars show standard errors of the mean. The average values over the gray shaded regions (4.0−4.5 nm) were considered as the dissociation free energy ΔGd. The average values over this region are shown in Table 3.

typically used in MSMs constructed from microsecond trajectories for folding/unfolding studies. However, the fraction of the total simulation time in each trial versus τ is approximately 104, consistent with the value suggested for enhanced sampling methods to achieve convergence of MSMs.73 Moreover, Suarez et al. conducted detailed analysis of Mean First Passage Time (MFPT) with non-Markovian estimators and found a reduction in the bias intrinsic to Markov MFPT estimation, even at the shortest lag times or simple discretization of the configuration in one-dimensional space.74 In addition, Zhang et al. used 0.5 ps for τ, which is much shorter than our value, to build transition matrices for MSMs from replica exchange simulation.75 Therefore, we believe our choice of τ is reasonable. Since microstates were assigned in 1D space, transition matrix elements Tij were relatively sparse, and only neighboring elements had significantly large values if index of microstates (i and j of Tij) are rank-ordered along the COM distance. In the case of PaCS-MD100,0.1, average of elements Tii, Tii±1, Tii±2, and Tii±3 were 0.33, 0.11, 0.05, and 0.03, which indicates that 71% of transitions occurred within the third neighbors along COM distance. The obtained dissociation free energies were 27.8 ± 0.8 and 30.5 ± 0.8 kJ/mol for PaCSMD10,1 and PaCS-MD100,0.1, respectively (Table 3), which are comparable to the binding free energy values measured by isothermal titration calorimetry (ITC) and surface plasmon resonance (SPR).32,76 412

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Article

Journal of Chemical Theory and Computation

The following analysis (see next section) shows significant differences in dissociation processes between PaCS-MD and SMD. The steering force in SMD induced a biased dissociation process, which might significantly change the dissociation free energy. 3.7. Disruption of LYZ-triNAG Interactions during the Dissociation Process. We examined the dissociation process of triNAG from LYZ by analyzing the breakage of the intermolecular hydrogen bonds shown in Figure 1(d,e) and determined the order in which LYZ residues dissociated from triNAG in each trial (Tables S1 and S2). There are 5.6 hydrogen bonds between LYZ and triNAG in the bound state during 1 μs conventional MD, and additional transient hydrogen bonds were formed in the partially bound state. The number of hydrogen bonds between LYZ and triNAG during dissociation in PaCS-MD10,0.1, PaCS-MD100,0.1, and PaCS-MD10,1 simulations were on average 8.1, 8.7, and 10.6, and those in SMDfast, SMDmed, and SMDslow simulations were 7.5, 9.4, and 12.3, respectively. These results indicate that the longer the MD run, or the slower the SMD velocity, the larger the number of hydrogen bonds formed during dissociation. We found that ASN59, TRP62, TRP63, and ALA107 are key residues that always formed hydrogen bonds with triNAG in the bound state. Interestingly, the TRP62 and TRP63 hydrogen bonds with NAG III located in pocket C broke before those of ASN59 and ALA107 with NAG II situated in pocket C, regardless of the simulation type. ASN59, TRP62, and TRP63 belong to the β domain, while ALA107 is situated on the α domain and is essentially outside the cleft (Figure 1(d,e)). NAG II, which binds to ASN103 and ASP101 of pocket B, dissociated after NAG III. triNAG did not slide along the binding pockets during dissociation but rather dissociated perpendicularly from the binding cleft starting from NAG III. If such sliding indeed happens, it should simultaneously break all the intermolecular hydrogen bonds, which should suddenly increase the free energy. The order of hydrogen bond breakage during SMD depended on the pulling velocity; however, no clear variations were seen during PaCS-MD. SMD can lead to the dissociation orders, being different from those observed in PaCS-MD. ASP48 is situated farther and deeper in the binding cleft compared to ASN59 and rarely forms hydrogen bonds with triNAG in the equilibrium bound state. However, during dissociation, ASP48 frequently formed hydrogen bonds with triNAG and then broke during SMDmed (9/12 trials) and SMDslow (11/12 trials). In PaCS-MD10,1 simulations, ARG112 formed transient hydrogen bonds with triNAG, and ASN106 played the same role in PaCS-MD100,0.1 simulation. ASN106 and ARG112 are located on the surface of the α domain. The hydrogen bond between ARG73 and NAG III tended to break first during PaCS-MD, whereas that between ARG61 and NAG I was the first hydrogen bond lost during SMD. These hydrogen bonds were both formed during dissociation. In PaCS-MD simulation, NAG I tended to dissociate first and was exposed to water before NAG III dissociation, whereas these events were reversed during SMD. In all the cases of SMD, NAG II is found to dissociate after both NAG I and NAG III do. The difference in dissociation order in SMD simulation might be due to force bias, which may contribute to a higher dissociation free energy in SMDfast/Jarzynski.

Table 3. Comparison of Dissociation and Association Free Energies Obtained by Simulations and Experiments quantities ΔGd

ΔGb

methods 10,1

PaCS-MD /MSM PaCS-MD100,0.1/MSM PaCS-MD10,1/US PaCS-MD10,0.1/US SMDfast/US SMDmed/US SMDslow/US SMDfast/Jarzynski SMDmed/Jarzynski SMDslow/Jarzynski ITC at pH 4.632 ITC at pH 7.332 SPR at pH 7.432 ITC at pH 7.076

free energy differencea (kJ/mol) 27.8 30.5 37.7 26.8 29.7 27.2 22.6 148.1 97.5 67.4

± ± ± ± ± ± ± ± ± ±

0.8 (27.6 ± 0.8) 0.8 (31.3 ± 0.8) 1.9 (39.8 ± 3.2) 1.3 (26.6 ± 1.2) 0.9 (30.1 ± 0.9) 1.1 (27.4 ± 1.1) 0.9 (23.0 ± 1.0) 2.3 (142.1 ± 2.1) 2.2 (98.0 ± 2.0) 1.5 (66.9 ± 1.5) −29.2 −28.5 −26.9 −28.9

a Average values over d = 4.0−4.5 nm. The values after ± show standard errors of the mean. The values in parentheses indicate averages and standard error at 4.0 nm only over trials. Note that the simulations were conducted at 300 K, and all the experiments were performed at 25 °C.

We also calculated the free energy profile using PaCS-MD/ US (Figure 7(b)). Sharp peaks observed at around 1.6 nm were not observed in the PaCS-MD/MSM results. While the free energy profile in PaCS-MD/MSM was directly calculated from PaCS-MD trajectories, that in PaCS-MD/US was obtained by additional multiple umbrella samplings with restrained interCOM distances. The peak at 1.6 nm indicates the region where triNAG was trapped in the binding pockets and strongly corresponds to the first peak of steering forces in SMD. The minimum at 1.9 nm after the first peak is related to the upper limit of the bound state, as shown in Figure 2. PaCS-MD10,1/ US yielded a free energy difference larger than the experimental values, while PaCS-MD10,0.1/US gave a more reasonable value, 26.8 ± 1.3 kJ/mol. The energy profiles of SMD/US indicated similar tendencies to those of PaCS-MD/US, including the position of the local minimum at 1.9 nm (Figure 7(c)). This minimum was clearly seen in SMDfast and SMDmed simulations but was not observed in SMDslow simulations. The binding free energies obtained by SMDfast/US (−29.7 ± 0.9 kJ/mol) and SMDmed/US (−27.2 ± 1.1 kJ/mol) are in the range of available experimental results (Table 3), whereas SMDslow/US underestimated ΔGb. SMD/Jarzynski (Figure 7(d)) significantly overestimated the dissociation free energy between LYZ and triNAG. Yamashita and Fujitani showed that SMD/US yielded higher free energy difference than did US with multistep targeted MD, which applied restraints to the heavy atoms of proteins with harmonic potential in the lysozyme-HyHEL complex.24 We did not observe such overestimation in SMD/US, probably because the conformational change of the small ligand can be restored during US. The dissociation free energy profiles obtained from SMD/Jarzynski showed a clear tendency that lower SMD pulling velocities resulted in lower binding free energy (Figure 7(d)). The first plateaus around 1.5−2.1 nm correspond to the region between the first and second force peaks shown in Table 2. The heights of the plateaus were 93.3 (SMDfast), 46.9 (SMDmed), and 20.5 kJ/mol (SMDslow). The dissociation free energies obtained by SMD/Jarzynski (Table 3) were larger by factors of 2−5 than that expected from experimental results.32

4. CONCLUSION AND DISCUSSION PaCS-MD is a straightforward simulation algorithm, as is its easy implementation, and is suitable for parallel and/or 413

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Article

Journal of Chemical Theory and Computation distributed computing. In this work, we first showed that PaCSMD easily induced the dissociation of LYZ and triNAG within the order of 100−101 ns and total cost of 102 ns MD simulation time (Table 1 and Figure 1). In contrast, no dissociation was observed during 1 μs of conventional MD (Figure S1). The cycles of multiple short MD simulations and the selection of rare events clearly accelerated dissociation. Although no bias was applied during MD, the dissociation speed of triNAG in PaCS-MD10,1 was equivalent to the pulling velocity of SMDfast, whereas those in PaCS-MD10,0.1 and PaCS-MD100,0.1 were 5 times faster. It is natural to ask whether PaCS-MD can be used to dissociate other ligands, especially larger molecules. Since the binding free energy of small ligands can be calculated using other methods, such as TI, we plan to apply the PaCS-MD/ MSM method to combinations of two larger molecules, which are much more challenging for all-atom binding free energy calculations. We consider the present work as the first step toward all-atom binding free energy calculations of protein− protein and other large molecule complexes. Thus, we focused only on LYZ-triNAG and examined many aspects of the dissociation process very carefully in this paper. We are now conducting dissociation PaCS-MD simulations of protein fragments and intact proteins from other proteins. To date, all dissociations were induced easily and will be reported in the future. We examined the effects of the number of replicas (nrep), MD length (tcyc), velocity rerandomization, and snapshot selection on PaCS-MD sampling. As clearly shown by comparison of the results obtained using PaCS-MD100,0.1 and PaCS-MD10,1 (Table 1 and Figure 2), increasing nrep is a more efficient way to reduce the number of cycles necessary for dissociation than increasing tcyc with a given fixed total simulation time per cycle, because the probability of observing rare events is proportional to nrep. We found that velocity rerandomization enhanced sampling toward dissociation in the bound state, in which triNAG was trapped in energy minima by interaction with LYZ. In this situation, velocity rerandomization acted as a perturbation to enhance the occurrence of fluctuations toward escape from the bound state and selection raised the probability of a rare event occurrence. Diffusion plays a more important role in the unbound state. Diffusion-governed movements occurring near the end of each MD run tend to be selected in PaCS-MD, which accelerates the dissociation process. One trial of PaCS-MD generates a dissociation pathway as a combination of multiple short MD trajectories that mutually overlap in conformational space. We confirmed that the generated trajectories can be properly utilized to construct a MSM. The PaCS-MD-generated pathway is not a true dissociation pathway as a function of time, yet it can provide statistical information regarding the dissociation process along a zigzag-rod-like conformational space that connects the bound state to the completely unbound state. Interestingly, a drastic (10-fold) increase of nrep in PaCS-MD100,0.1 compared to PaCSMD10,0.1 only resulted in a slight increase in σ each (d)(Figure 6), indicating that PaCS-MD100,0.1 sampled a conformational space similar to that in PaCS-MD10,0.1 but more densely. This is related to the fact that we could build a MSM from each trial of PaCS-MD100,0.1, whereas the statistics were insufficient in PaCSMD10,0.1. The longer (10-fold) simulation time, tcyc, in PaCSMD10,1 generated an obviously wider conformational space along the dissociation pathway per trial, as shown by the value

of σ each (d) more than doubling. Since PaCS-MD100,0.1 and PaCS-MD10,1 required the same computational time, PaCSMD100,0.1 densely sampled a narrower conformational space along the pathway, while PaCS-MD10,1 explored a wider space more sparsely. PaCS-MD10,1/MSM resulted in a dissociation free energy ΔGd slightly closer to the experimental value, but the effect was relatively small. As long as sufficient statistics are achieved by sampling along the pathway, the width of the pathway, which shows the range of sampling almost perpendicular to the pathway, should not significantly affect ΔGd. Although ΔGd is in principle independent of the pathway, the actual value obtained by each trial likely contains some errors. Therefore, ΔGd was calculated as the average over multiple pathways. Note that the average in the Jarzynski equality is taken over much narrower dissociation pathways generated by SMD. In this case, greater statistical accuracy can be achieved only by obtaining more trajectories. To obtain a more accurate ΔGd, averaging over multiple trials may be more important than generating a wider dissociation pathway. To build a MSM directly from PaCS-MD trajectories, we should carefully choose nrep and tcyc. As discussed above, nrep affects both statistics and the number of PaCS-MD cycles required to achieve dissociation. The number of cycles can be reduced more efficiently by increasing nrep. Building a proper MSM also requires the correct choice of tcyc: tcyc affects the choice of lag time τ and the definition of microstates in MSM, because the condition τ ≤ tcyc/2 is expected from the statistical point of view.45 Therefore, the merit of using longer τ is to provide flexibility in deciding proper τ and microstates. In this study, as mentioned above, we judged that tcyc = 0.1 ns in PaCS-MD100,0.1 simulations was sufficient to build a MSM. It should be noted that we obtained the comparable ΔGd values from PaCS-MD100,0.1 and PaCS-MD10,1 simulations despite using different lag times. Using PaCS-MD to generate initial pathways with a limited number of nrep, followed by MD simulations to obtain more statistics to build a MSM, provides more options. In this case, PaCS-MD can be mainly dedicated to accelerating expected movements with sufficiently short MDs. A MSM can be built by extending the MD simulation length and/or by adding more MD simulations later. Also, as in the case where we noticed insufficient statistics to build a MSM after PaCS-MD, more MD trajectories can be added to properly construct a MSM. As noted above, each trial may not provide adequate statistics for MSM. One approach to addressing this might be to gather trajectories from multiple trials and build one MSM. We constructed MSMs for all trajectories of the ten trials conducted using PaCS-MD10,0.1 and obtained ΔG = 27.4(kJ/ mol), consistent with experimental data. However, it should be noted that snapshots with similar d values can be assigned to one microstate, but snapshots from different trials can be very far in 3D space, especially in the unbound state, and thus should not be categorized in the same microstate. This is clearly seen in Figure S3(a). Red small spheres in the unbound state in Figure S3(a) are very far from green spheres in the unbound state in 3D space; however, they can be categorized into the same microstate if d is very close. We judged that using d values to define microstates is invalid from the physical point of view. To overcome this problem, we would suggest using 3D COM coordinates to build MSM and convert the result to PMF along 1D reaction coordinate PMF (inter-COM distance), rather than building MSM based on the 1D data set. 414

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Journal of Chemical Theory and Computation As shown in Table 3, ΔGd calculation with PaCS-MD/MSM yielded better results than calculation with PaCS-MD/US, SMD/US, or SMD/Jarzynski, as mentioned above. In addition, two different PaCS-MD/MSMs gave similar ΔGd results. PaCSMD/MSM showed the additional advantage that it gave the lowest standard deviation, which indicates the least free energy variation among the trials. Notably, the total simulation time for PaCS-MD/MSM is less than PaCS-MD/US or SMD/US, because US requires longer simulation time (Table 1). The ΔGd values from SMD/US and SMD/Jarzynski clearly depended on the SMD pulling velocity. In particular, SMD/ Jarzynski overestimated ΔGd in all cases. Detailed analysis of the LYZ-triNAG interactions in SMD showed that the dissociation order of the three NAGs from LYZ was different from that in PaCS-MD, suggesting that an unnatural dissociation process was induced by the force bias in SMD. The results of SMD/US indicated that artifacts caused by this bias were mostly recovered during US, but velocity dependence was still observed. Therefore, the choice of pulling velocity is important in SMD/US simulations.

■

■

ACKNOWLEDGMENTS

■

REFERENCES

D.T. and A.K. thank Jacob Swadling for carefully reading the manuscript.

(1) Lelièvre, T.; Rousset, M.; Stoltz, G. Free Energy Computations; Imperial College Press: London, 2010. (2) Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. J. Chem. Phys. 1935, 3 (5), 300−313. (3) Zwanzig, R. Nonlinear Generalized Langevin Equations. J. Stat. Phys. 1973, 9 (3), 215−220. (4) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. THE Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. J. Comput. Chem. 1992, 13 (8), 1011−1021. (5) Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23 (2), 187−199. (6) Jarzynski, C. Nonequilibrium Equality for Free Energy Differences. Phys. Rev. Lett. 1997, 78 (14), 2690−2693. (7) Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2003, 90 (23), 238302. (8) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99 (20), 12562−12566. (9) Wang, F.; Landau, D. P. Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States. Phys. Rev. Lett. 2001, 86 (10), 2050−2053. (10) Jorgensen, W.; Ravimohan, C. Monte Carlo Simulation of Differences in Free Energies of Hydration. J. Chem. Phys. 1985, 83 (6), 3050−3054. (11) Straatsma, T. P.; McCammon, J. A. Multiconfiguration Thermodynamic Integration. J. Chem. Phys. 1991, 95 (2), 1175. (12) Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. Free Energy Calculation from Steered Molecular Dynamics Simulations Using Jarzynski’s Equality. J. Chem. Phys. 2003, 119 (6), 3559−3566. (13) Suan Li, M.; Khanh Mai, B. Steered Molecular Dynamics-A Promising Tool for Drug Design. Curr. Bioinf. 2012, 7 (4), 342−351. (14) Ramírez, C. L.; Martí, M. A.; Roitberg, A. E. Chapter Six − Steered Molecular Dynamics Methods Applied to Enzyme Mechanism and Energetics. Methods Enzymol. 2016, 578, 123−143. (15) Sugita, Y.; Okamoto, Y. Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 1999, 314 (1), 141− 151. (16) Bernardi, R. C.; Melo, M. C. R.; Schulten, K. Enhanced Sampling Techniques in Molecular Dynamics Simulations of Biological Systems. Biochim. Biophys. Acta, Gen. Subj. 2015, 1850 (5), 872−877. (17) Harada, R.; Kitao, A. Parallel Cascade Selection Molecular Dynamics (PaCS-MD) to Generate Conformational Transition Pathway. J. Chem. Phys. 2013, 139 (3), 035103. (18) Deng, Y.; Roux, B. Computations of Standard Binding Free Energies with Molecular Dynamics Simulations. J. Phys. Chem. B 2009, 113 (8), 2234−2246. (19) Gumbart, J. C.; Roux, B.; Chipot, C. Standard Binding Free Energies from Computer Simulations: What Is the Best Strategy? J. Chem. Theory Comput. 2013, 9 (1), 794−802. (20) Deng, Y.; Roux, B. Calculation of Standard Binding Free Energies: Aromatic Molecules in the T4 Lysozyme L99A Mutant. J. Chem. Theory Comput. 2006, 2 (5), 1255−1273. (21) Gumbart, J. C.; Roux, B.; Chipot, C. Efficient Determination of Protein−Protein Standard Binding Free Energies from First Principles. J. Chem. Theory Comput. 2013, 9 (8), 3789−3798. (22) Sugita, Y.; Kitao, A.; Okamoto, Y. Multidimensional ReplicaExchange Method for Free-Energy Calculations. J. Chem. Phys. 2000, 113 (15), 6042−6051.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00504. Residual order of hydrogen bond breakage with triNAG in PaCS-MD and SMD (Tables S1 and S2); root-meansquared-deviation of LYZ, triNAG, and inter-COM distance between them during 1 μs conventional MD simulation (Figure S1); mean-squared-displacement of triNAG calculated for the reactive trajectories of PaCSMD (Figure S2); visualization of the center-of-mass position of triNAG along dissociation pathways (Figure S3) (PDF)

■

Article

AUTHOR INFORMATION

Corresponding Author

*Phone: +81-3-5734-3373. Fax: +81-3-5734-3372. E-mail: [email protected]. ORCID

Akio Kitao: 0000-0002-5221-0806 Funding

This research was supported by MEXT/JSPS KAKENHI (Nos. 25104002 and 15H04357) to A.K. and by MEXT as “Priority Issue on Post-K Computer” (Building Innovative Drug Discovery Infrastructure Through Functional Control of Biomolecular Systems) to A.K. The computations were partly performed using the supercomputers at the RCCS, The National Institute of Natural Science, and ISSP, The University of Tokyo. This research also used computational resources of the K computer and other computers of the HPCI system provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Project ID: hp140031, hp150049, hp150270, hp160207, and hp170254). D.T. acknowledges financial support from Hitachi Scholarship program No. S175 of the Hitachi Global Foundation. Notes

The authors declare no competing financial interest. 415

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Article

Journal of Chemical Theory and Computation (23) Schlitter, J.; Engels, M.; Krüger, P. Targeted Molecular Dynamics: A New Approach for Searching Pathways of Conformational Transitions. J. Mol. Graphics 1994, 12 (2), 84−89. (24) Yamashita, T.; Fujitani, H. On Accurate Calculation of the Potential of Mean Force between Antigen and Antibody: A Case of the HyHEL-10-Hen Egg White Lysozyme System. Chem. Phys. Lett. 2014, 609, 50−53. (25) Harada, R.; Kitao, A. Nontargeted Parallel Cascade Selection Molecular Dynamics for Enhancing the Conformational Sampling of Proteins. J. Chem. Theory Comput. 2015, 11 (11), 5493−5502. (26) Harada, R.; Nakamura, T.; Shigeta, Y. Sparsity-Weighted Outlier FLOODing (OFLOOD) Method: Efficient Rare Event Sampling Method Using Sparsity of Distribution. J. Comput. Chem. 2016, 37 (8), 724−738. (27) Jani, V.; Sonavane, U.; Joshi, R. Traversing the Folding Pathway of Proteins Using Temperature-Aided Cascade Molecular Dynamics with Conformation-Dependent Charges. Eur. Biophys. J. 2016, 45 (5), 463−482. (28) Peng, J.; Zhang, Z. Unraveling Low-Resolution Structural Data of Large Biomolecules by Constructing Atomic Models with Experiment-Targeted Parallel Cascade Selection Simulations. Sci. Rep. 2016, 6 (1), 29360. (29) Carrillo, W.; Garcia-Ruiz, A.; Recio, I.; Moreno-Arribas, M. V. Antibacterial Activity of Hen Egg White Lysozyme Modified by Heat and Enzymatic Treatments against Oenological Lactic Acid Bacteria and Acetic Acid Bacteria. J. Food Prot. 2014, 77 (10), 1732−1739. (30) Ming, D.; Wall, M. E. Quantifying Allosteric Effects in Proteins. Proteins: Struct., Funct., Genet. 2005, 59 (4), 697−707. (31) Cheetham, J. C.; Artymiuk, P. J.; Phillips, D. C. Refinement of an Enzyme Complex with Inhibitor Bound at Partial Occupancy: Hen Egg-White Lysozyme and Tri-N-Acetylchitotriose at 1.75 Å Resolution. J. Mol. Biol. 1992, 224 (3), 613−628. (32) Ishikawa, T.; Burri, R. R.; Kamatari, Y. O.; Sakuraba, S.; Matubayasi, N.; Kitao, A.; Kuwata, K. A Theoretical Study of the Two Binding Modes between Lysozyme and Tri-NAG with an Explicit Solvent Model Based on the Fragment Molecular Orbital Method. Phys. Chem. Chem. Phys. 2013, 15 (10), 3646. (33) Takemura, K.; Burri, R. R.; Ishikawa, T.; Ishikura, T.; Sakuraba, S.; Matubayasi, N.; Kuwata, K.; Kitao, A. Free-Energy Analysis of lysozyme−triNAG Binding Modes with All-Atom Molecular Dynamics Simulation Combined with the Solution Theory in the Energy Representation. Chem. Phys. Lett. 2013, 559, 94−98. (34) Zhong, Y.; Patel, S. Binding Structures of Tri- N -Acetyl-βGlucosamine in Hen Egg White Lysozyme Using Molecular Dynamics with a Polarizable Force Field. J. Comput. Chem. 2013, 34 (3), 163− 174. (35) Izrailev, S.; Stepaniants, S.; Isralewitz, B.; Kosztin, D.; Lu, H.; Molnar, F.; Wriggers, W.; Schulten, K. Steered Molecular Dynamics. In Computational Molecular Dynamics: Challenges, Methods, Ideas SE 2; 1999; Vol. 4, pp 39−65. (36) Sotomayor, M.; Schulten, K. Single-Molecule Experiments in Vitro and in Silico. Science (Washington, DC, U. S.) 2007, 316 (5828), 1144. (37) 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. (38) Chodera, J. D.; Noé, F. Markov State Models of Biomolecular Conformational Dynamics. Curr. Opin. Struct. Biol. 2014, 25, 135−144. (39) Kitao, A.; Hirata, F.; Go̅, N. The Effects of Solvent on the Conformation and the Collective Motions of Protein: Normal Mode Analysis and Molecular Dynamics Simulations of Melittin in Water and in Vacuum. Chem. Phys. 1991, 158 (2−3), 447−472. (40) Molgedey, L.; Schuster, H. G. Separation of a Mixture of Independent Signals Using Time Delayed Correlations. Phys. Rev. Lett. 1994, 72 (23), 3634−3637. (41) Pérez-Hernández, G.; Paul, F.; Giorgino, T.; De Fabritiis, G.; Noé, F. Identification of Slow Molecular Order Parameters for Markov Model Construction. J. Chem. Phys. 2013, 139 (1), 015102.

(42) Naritomi, Y.; Fuchigami, S. Slow Dynamics in Protein Fluctuations Revealed by Time-Structure Based Independent Component Analysis: The Case of Domain Motions. J. Chem. Phys. 2011, 134 (6), 065101. (43) 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. (44) Prinz, J. H.; Wu, H.; Sarich, M.; Keller, B.; Senne, M.; Held, M.; Chodera, J. D.; Schütte, C.; Noé, F. Markov Models of Molecular Kinetics: Generation and Validation. J. Chem. Phys. 2011, 134 (17), 174105. (45) Doerr, S.; de Fabritiis, G. On-the-Fly Learning and Sampling of Ligand Binding by High-Throughput Molecular Simulations. J. Chem. Theory Comput. 2014, 10 (5), 2064−2069. (46) Beauchamp, K. A.; Bowman, G. R.; Lane, T. J.; Maibaum, L.; Haque, I. S.; Pande, V. S. MSMBuilder2: Modeling Conformational Dynamics on the Picosecond to Millisecond Scale. J. Chem. Theory Comput. 2011, 7 (10), 3412−3419. (47) Grossfield, A. WHAM; Grossfield Lab. http://membrane.urmc. rochester.edu/content/wham (accessed Jan 4, 2017). (48) Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. Free Energy Calculation from Steered Molecular Dynamics Simulations Using Jarzynski’s Equality. J. Chem. Phys. 2003, 119 (6), 3559−3566. (49) Dellago, C.; Hummer, G. Computing Equilibrium Free Energies Using Non-Equilibrium Molecular Dynamics. Entropy 2014, 16 (1), 41−61. (50) Park, S.; Schulten, K. Calculating Potentials of Mean Force from Steered Molecular Dynamics Simulations. J. Chem. Phys. 2004, 120 (13), 5946−5961. (51) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Struct., Funct., Genet. 2006, 65, 712−725. (52) Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; GonzalezOuteirino, J.; Daniels, C. R.; Foley, B. L.; Woods, R. J. GLYCAM06: A Generalizable Biomolecular Force Field. Carbohydrates. J. Comput. Chem. 2008, 29 (4), 622−655. (53) Pronk, S.; Pall, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29 (7), 845−854. (54) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18 (12), 1463−1472. (55) Van Gunsteren, W. F.; Berendsen, H. J. C. A Leap-Frog Algorithm for Stochastic Dynamics. Mol. Simul. 1988, 1 (3), 173−185. (56) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. A Computer Simulation Method for the Calculation of Equilibrium Constants for the Formation of Physical Clusters of Molecules: Application to Small Water Clusters. J. Chem. Phys. 1982, 76 (1), 637− 649. (57) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126 (1), 014101. (58) Nosé, S. A. Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81 (1), 511−519. (59) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31 (3), 1695− 1697. (60) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, a.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81 (1984), 3684−3690. (61) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52 (12), 7182−7190. (62) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit Reversible Integrators for Extended Systems Dynamics. Mol. Phys. 1996, 87 (5), 1117−1157. 416

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417

Article

Journal of Chemical Theory and Computation (63) Yeh, I.-C.; Hummer, G. Diffusion and Electrophoretic Mobility of Single-Stranded RNA from Molecular Dynamics Simulations. Biophys. J. 2004, 86 (2), 681−689. (64) Yeh, I. C.; Hummer, G. System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions. J. Phys. Chem. B 2004, 108 (40), 15873−15879. (65) Takemura, K.; Kitao, A. Effects of Water Model and Simulation Box Size on Protein Diffusional Motions. J. Phys. Chem. B 2007, 111 (41), 11870−11872. (66) Matsumura, I.; Kirsch, J. F. Is Aspartate 52 Essential for Catalysis by Chicken Egg White Lysozyme? The Role of Natural Substrate-Assisted Hydrolysis. Biochemistry 1996, 35 (6), 1881−1889. (67) Noguchi, S.; Miyawaki, K.; Satow, Y. Succinimide and Isoaspartate Residues in the Crystal Structures of Hen Egg-White Lysozyme Complexed with Tri-N-Acetylchitotriose. J. Mol. Biol. 1998, 278 (1), 231−238. (68) Ohmura, T.; Ueda, T.; Ootsuka, K.; Saito, M.; Imoto, T. Stabilization of Hen Egg White Lysozyme by a Cavity-Filling Mutation. Protein Sci. 2001, 10 (2), 313−320. (69) Kamiya, N.; Yonezawa, Y.; Nakamura, H.; Higo, J. ProteinInhibitor Flexible Docking by a Multicanonical Sampling: Native Complex Structure with the Lowest Free Energy and a Free-Energy Barrier Distinguishing the Native Complex from the Others. Proteins: Struct., Funct., Genet. 2008, 70 (1), 41−53. (70) Patriksson, A.; van der Spoel, D. A Temperature Predictor for Parallel Tempering Simulations. Phys. Chem. Chem. Phys. 2008, 10 (15), 2073. (71) Saxton, M. J. Single-Particle Tracking: The Distribution of Diffusion Coefficients. Biophys. J. 1997, 72 (4), 1744−1753. (72) Horton, M.; Charras, G.; Lehenkari, P. ANALYSIS OF LIGAND−RECEPTOR INTERACTIONS IN CELLS BY ATOMIC FORCE MICROSCOPY. J. Recept. Signal Transduction Res. 2002, 22 (1−4), 169−190. (73) 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, 99−105. (74) Suárez, E.; Adelman, J. L.; Zuckerman, D. M. Accurate Estimation of Protein Folding and Unfolding Times: Beyond Markov State Models. J. Chem. Theory Comput. 2016, 12 (8), 3473−3481. (75) Zhang, B. W.; Dai, W.; Gallicchio, E.; He, P.; Xia, J.; Tan, Z.; Levy, R. M. Simulating Replica Exchange: Markov State Models, Proposal Schemes, and the Infinite Swapping Limit. J. Phys. Chem. B 2016, 120 (33), 8289−8301. (76) Ogata, M.; Umemoto, N.; Ohnuma, T.; Numata, T.; Suzuki, A.; Usui, T.; Fukamizo, T. A Novel Transition-State Analogue for Lysozyme, 4-O-β-Tri-N-Acetylchitotriosyl Moranoline, Provided Evidence Supporting the Covalent Glycosyl-Enzyme Intermediate. J. Biol. Chem. 2013, 288 (9), 6072−6082. (77) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14 (1), 33−38. (78) Laskowski, R. A.; Swindells, M. B. LigPlot+: Multiple Ligand− Protein Interaction Diagrams for Drug Discovery. J. Chem. Inf. Model. 2011, 51 (10), 2778−2786.

417

DOI: 10.1021/acs.jctc.7b00504 J. Chem. Theory Comput. 2018, 14, 404−417