Integrated Hamiltonian Sampling: A Simple and Versatile Method for

Article pubs.acs.org/JPCB

Integrated Hamiltonian Sampling: A Simple and Versatile Method for Free Energy Simulations and Conformational Sampling Toshifumi Mori,*,†,∥ Robert J. Hamers,‡ Joel A. Pedersen,§ and Qiang Cui*,† †

Department of Chemistry and Theoretical Chemistry Institute, University of WisconsinMadison, 1101 University Avenue, Madison, Wisconsin 53706, United States ‡ Department of Chemistry, University of WisconsinMadison, 1101 University Avenue, Madison, Wisconsin 53706, United States § Departments of Soil Science, Civil & Environmental Engineering, and Chemistry, University of WisconsinMadison, 1525 Observatory Drive, Madison, Wisconsin 53706, United States S Supporting Information *

ABSTRACT: Motivated by specific applications and the recent work of Gao and co-workers on integrated tempering sampling (ITS), we have developed a novel sampling approach referred to as integrated Hamiltonian sampling (IHS). IHS is straightforward to implement and complementary to existing methods for free energy simulation and enhanced configurational sampling. The method carries out sampling using an effective Hamiltonian constructed by integrating the Boltzmann distributions of a series of Hamiltonians. By judiciously selecting the weights of the different Hamiltonians, one achieves rapid transitions among the energy landscapes that underlie different Hamiltonians and therefore an efficient sampling of important regions of the conformational space. Along this line, IHS shares similar motivations as the enveloping distribution sampling (EDS) approach of van Gunsteren and coworkers, although the ways that distributions of different Hamiltonians are integrated are rather different in IHS and EDS. Specifically, we report efficient ways for determining the weights using a combination of histogram flattening and weighted histogram analysis approaches, which make it straightforward to include many end-state and intermediate Hamiltonians in IHS so as to enhance its flexibility. Using several relatively simple condensed phase examples, we illustrate the implementation and application of IHS as well as potential developments for the near future. The relation of IHS to several related sampling methods such as Hamiltonian replica exchange molecular dynamics and λ-dynamics is also briefly discussed.

■

INTRODUCTION Carrying out adequate conformational sampling is of great importance to meaningful condensed phase simulations, especially for biomolecular systems.1−6 Following decades of efforts, many powerful sampling methodologies are available in the literature: some are motivated by locating low free energy basins,7−23 while others target free energy differences between different states.24−31 Depending on the problem at hand, and also the computing resources available, specific sampling methods are likely effective. For example, when many processors are available, multidimensional Hamiltonian replica exchange molecular dynamics (HREMD)14,15 and related variants29,32,33 can be very powerful in attacking complex ligand binding free energy calculations.34,35 Similarly, great success has been reported for small protein folding simulations using simulated tempering,16,17 replica exchange,18 and related methods.36 Despite these advances, new sampling/free energy methods are beneficial to specific types of applications. The development described here, for example, is motivated by our interest in two types of problems. In the first, we are interested in computing © 2014 American Chemical Society

the free energy difference between two levels of potential functions; in particular, one potential function is ab initio quantum mechanical/molecular mechanical (QM/MM37−41), and the other is given by either a low-level (e.g., semiempirical42−45) QM/MM potential or a pure MM force field.46 As shown in the thermodynamic cycle in Figure 1a, by computing the free energy difference between two levels of potential functions (i.e., the vertical paths, ΔGh→l A/B), one is able to obtain the free energy change for a process with a high-level l h→l + ΔGh→l QM/MM potential (i.e., ΔGhA→B = ΔGA→B B − ΔGA ). It is advantageous to follow the thermodynamic cycle in Figure 1 because the most extensive sampling is carried out at low level l ), since the vertical paths (ΔGh→l (ΔGA→B A/B) implicate smaller changes and therefore should converge fairly quickly. In reality, however, when a high-level QM/MM method is used, even Special Issue: James L. Skinner Festschrift Received: February 7, 2014 Revised: March 11, 2014 Published: March 18, 2014 8210

dx.doi.org/10.1021/jp501339t | J. Phys. Chem. B 2014, 118, 8210−8220

The Journal of Physical Chemistry B

Article

pioneered by van Gunsteren and co-workers;63−65 in both approaches, one seeks an effective Hamiltonian66 that has good overlapping distributions with two “end states” and efficient sampling is carried out using this effective Hamiltonian. The ways that such an effective Hamiltonian is constructed, however, are rather different in IHS and EDS, and we believe IHS offers potentially a more general and flexible solution. Finally, IHS is also closely related to HREMD,14,15 although instead of having many parallel replicas each with a different Hamiltonian, IHS runs a single simulation using an effective Hamiltonian; the comparison between the two and with other related methods such as λ-dynamics24 is further elaborated below. In the following, we first describe the key idea and implementation of IHS. We then use several relatively simple examples (see Figure 2) to illustrate the effectiveness and Figure 1. Thermodynamic cycles used to compute (a) the free energy change of a process (A → B) at different (high, low) levels of theory and (b) the binding affinity of a molecule (e.g., a peptide) to a solid surface in solution. In each case, the quantity of interest is highlighted by the shaded notation, while the actual simulations are carried out along alternative thermodynamic pathways for higher computational efficiency. Quantities related to examples studied in the current work are colored in blue.

simulation of the 100 ps time scale is very demanding; thus, it is still worthwhile to develop sampling methods that minimize the computational cost associated with the computation of ΔGh→l A/B. Along this line, we note the recent work of Woodcock and coworkers,47 who demonstrated the value of using NonBoltzmann Bennett reweighing schemes for computing precisely the quantities of interest in Figure 1a. When the potential functions associated with the different levels do not overlap well, however, it is likely that including sampling at both levels is required. Another application that motivated our development is the study of peptide/protein binding at solid/liquid interfaces,48−58 a general and important problem prevalent to both biology and nano technology. As touched upon in our recent work,59 which discussed different strategies for computing the binding affinity of molecules to a solid surface (for example, Figure 1b), properly sampling the conformation of even a short peptide at the surface can be challenging because potentially strong local peptide/surface interaction can easily trap the system in local conformational minima. Therefore, one needs to enhance the sampling of both intrapeptide and peptide-surface degrees of freedom; the basic interfacial structure, however, should not be perturbed tremendously by the enhanced sampling protocol,60 suggesting that manipulating specific components of the Hamiltonian (or potential function) rather than global properties such as temperature is more appropriate in this context; alternatively, the temperature of specific degrees of freedom can be elevated by adiabatically decoupling them with the rest of the degrees of freedom.19,30 As explained in detail below, our specific development is largely inspired by the approach of integrated tempering sampling (ITS) of Gao.61,62 Rather than integrating distributions from different temperatures, we integrate distributions from different Hamiltonians (potential functions) to achieve efficient sampling; thus, we term our approach integrated Hamiltonian sampling (IHS) to honor the inspiration from ITS. We also recognize that our method shares key elements with the enveloping distribution sampling (EDS) approach

Figure 2. Several simple condensed phase examples used to illustrate the implementation and application of IHS: a dipole inversion process in solution; relative solvation free energy of a QM and a MM water in bulk water; conformational sampling of a short RGD peptide in solution and at the (110) surface of rutile TiO2 surface.

versatility of IHS. Finally, we make a few concluding remarks, especially concerning the comparison of IHS to several related methodologies such as HREMD.

■

THEORETICAL METHOD Integrated Hamiltonian Sampling (IHS). We set up our problem by defining two “end-state” Hamiltonians with potential functions U0 and U1, respectively. For free energy simulations, these represent the two physical states of interest; alternatively, U1 may represent an unphysical Hamiltonian (e.g., with torsional barriers flattened) introduced to facilitate the sampling of the physical system described by U0. As schematically illustrated in Figure 3, the two potential functions in general have global minima at different locations and therefore barriers need to be overcome to ensure that important free energy minima are sampled in both states. In general, the situation is complex because the nature of the essential coordinates along which significant barriers exist for the transition between U0 and U1 is not known. Overcoming such “hidden barriers” poses an important challenge to (bio)molecular simulations in general and to free energy simulations in particular. To facilitate transitions between U0 and U1, we introduce a series of intermediate states with potentials Uλ described by 8211



Article

the end states can be computed from IHS in a straightforward manner. Although our work is motivated by ITS, we note that the idea of sampling the configuration space with an effective Hamiltonian, defined in terms of the sum of distributions of two or more end states, was first proposed by C. Oostenbrink and W. F. van Gunsteren66 and was further developed as the envelope distribution sampling (EDS) method.64,65 In the EDS method, only the end states (U0 and U1) with weights are used, and in addition, a smoothness parameter was introduced to enhance the sampling. Although EDS has only a small number of parameters, determining their values was not easy, since the smoothness parameter and weight functions are not independent; the small number of parameters also may limit the effectiveness of EDS for complex systems, although a range of interesting applications have already been reported.63−65,67−69 By comparison, the weights are independent of each other (apart from a constant) in the IHS framework. This can be advantageous in that instead of running a single long trajectory we can easily combine samples from multiple short trajectories using the weighted histogram analysis method (WHAM)70 or the multistate Benett acceptance ratio (MBAR)71 with different weights, and the ease of using a large number of intermediates offers considerable flexibility to IHS. Moreover, generalizing IHS to multiple end states is straightforward without additional complication to the weight determination protocols. Below we first discuss how to update the weight function Ωi efficiently on the f ly and then the extension of IHS to multiple end state cases. Weight Optimization Protocols. The aim is to weight the contributions from different intermediate Hamiltonians such that transition among them, including the end states, is facile. To achieve this, one criterion is to set the expectation value of the probability for the weighted intermediate states to be uniform, or equivalently,

Figure 3. Schematic illustration of the IHS method. Blue and green lines represent the U0 and U1 “end-state” potentials, respectively, black and red lines represent the integrated potential, Ub (eq 3), without and with the weight functions (Ωi), and the thin lines between U0 and U1 are the intermediate potentials Uλi (eq 1). Ub and Ub′ are shifted in the −y direction to increase clarity. By judiciously choosing the weight functions (Ωi), important regions of both end states can be adequately sampled.

Uλ(R) = (1 − λ)U0(R) + λU1(R)

(1)

In standard free energy simulations, including λ dynamics, or HREMD simulations,14,15 these intermediate states are sampled separately. In IHS, we introduce an effective potential (Ub(R)) whose canonical distribution is the integrated distributions of the intermediate and end states 24

Ub(R) = −

1 ln β

∫0

1

dλ Ω(λ)e−βUλ(R)

(2)

Here Ω(λ) is a weight function to be determined. As schematically illustrated in Figure 3, a judicious choice of the weight function can lead to enhanced sampling of regions important to both end states. The details of the weight optimization procedure are described in the next section. In practice, the integral in eq 2 is discretized into a sum over Nλ windows N

Ub(R) = −

λ 1 ln ∑ Ωie−βUi(R) β i=1

−β −1 ln⟨Ωie−β(Ui − Ub)⟩b = −β −1 ln⟨Ωje−β(Uj − Ub)⟩b

Here Ui = Uλi, Ωi = Ωλi, and Nλ is the number of grids in the λ space. Note that, in the HREMD approach, where each replica computes a distinct potential function Ui, the computational cost increases linearly as a function of Nλ. In IHS, however, the number of intermediates can be increased with negligible overhead, especially when the intermediate potential is described by a linear combination of the end states (as in eq 1). Thus, the IHS approach allows one to use a very large number of intermediate states to optimize convergence. This is particularly advantageous when expensive potential functions such as high-level QM/MM potentials are used. From the probability distribution ρb(R) obtained via configurational sampling under Ub, the probability distribution of the weighted and original states i, ρ(b) i (R) and ρi(R), respectively, can be recovered via a reweighting procedure 1

ρi(b)(R) Zi(b)

=

1 ρ (R)Ωie−β(Ui(R) − Ub(R)) Zb b

1 1 ρ (R) = ρ (R)e−β(Ui(R) − Ub(R)) Zi i Zb b

(6)

Here ⟨···⟩b is the ensemble average over the configurations sampled with the potential Ub (eq 2). If we rewrite the weight function as Ũ i = β−1 ln Ωi, eq 6 can be rewritten as

(3)

Uĩ − Uj̃ = −β −1 ln

⟨e−β(Ui − Ub)⟩b ⟨e−β(Uj − Ub)⟩b

(7)

Uĩ = −β −1 ln⟨e−β(Ui − Ub)⟩b + C0 = −β −1 ln Zi + C0

(8)

or

with Zi and C0 being the (scaled) configuration integral of state i and an arbitrary constant, respectively. In the current work, C0 is determined to fulfill ∑i Ωi = ∑i exp(βŨ i) = 1. From eq 8, we realize that the most efficient sampling is achieved when Ũ i is the free energy of the unweighted state i. As mentioned above, the EDS approach determines the weights (and the smoothness parameters) from the information of a single simulation; in practice, this would require a long simulation (e.g., at least on the order of ∼ns). To make the IHS approach practical, it is desirable to estimate the weights quickly and also have the ability to refine them on the f ly, making use of all available information. To that end, here we explore the histogram flattening72 and WHAM70 approaches, and furthermore show that the combination of the two (referred to as “hybrid” below) can be very effective.

(4)

(5)

(b) Here Z(b) i = ∫ ρi dR and Zi = ∫ ρi dR are normalization factors (configurational integrals). Thus, any equilibrium properties of

8212



Article

In the histogram flattening approach, the weight for step n + 1 can be determined from the n-th sample via (n + 1) (n) (n) Uĩ = Uĩ + ΔUĩ

Ũ (n) i

Here C(n+1) is a normalization constant that can be omitted as we impose the constraint ∑i exp(βŨ (n+1) ) = 1. We see below i that WHAM is capable of updating the weights smoothly. Multistate IHS Method. Extension of IHS to include multiple end states is straightforward, yet has the potential advantage of enhancing multiple degrees of freedom simultaneously. In a similar manner to the multidimensional HREMD approaches,14 this can be done by increasing the dimensionality of λ, i.e., defining the intermediate states and optimizing the weight functions in the multidimensional λ space. In the case of three end states (U0, U1, and U2), for example, this can be given by

(9)

Ũ (n+1) i

with and being the weight for state i at step n and n + 1, respectively, and (n)

ΔUĩ

̃ (n)− Ub(n)) (n) ⟩b

= −β −1 ln⟨e−β(Ui − Ui

̃ (n)− Ub(n)) (n) ⟩b }j ]

− min[{−β −1 ln⟨e−β(Uj − Uj

(10)

While a similar approach has been used in ITS, since changes in the potential functions tend to be more drastic than temperature changes, a naı̈ve application of eq 10 to update the weights in IHS will lead to large oscillation in the estimated weights (e.g., see Figure 4). To avoid this issue and gently converge the weights, we apply a cap ΔUcap on the amount to be updated per step: (n + 1) (n) (n) Uĩ = Uĩ + min[ΔUĩ , ΔUcap]

Ub2D(R) = −

1 ln ∑ ∑ Ω λ1λ2 β λ1 λ 2

exp[−β {U0(R) + λ1(U1(R) − U0(R)) + λ 2(U2(R) − U0(R))}]

(14)

Here λ1 and λ2 are the indices for the two-dimensional (2D) λ space and Ωλ1λ2 is the weight in the 2D space to be determined. Alternatively, here we define the intermediate states as all the pairs of end states; e.g., in the three-end-state case, the integrated potential is described as

(11)

2

Ub(R) = − +

1 ln[ ∑ Ωme−βUm(R) + β m=0

∑ λi ≠ 0,1

Ωi e

−βUi12(R)

+

∑

01

Ωie−βUi

(R)

λi ≠ 0,1

∑ λi ≠ 0,1

20

Ωie−βUi

(R)

] (15)

where Umn i (R) = (1 − λi)Um(R) + λiUn(R). This has the advantage in that the λ space to be sampled is much less compared to the 2D case, and is expected to be even more beneficial when more end states are introduced.

■

EXAMPLES AND DISCUSSIONS Simulation Setup. As summarized in Figure 2, several relatively simple condensed phase problems are chosen to illustrate the implementation and application of IHS. In all calculations, the TIP3P water model73 and the AMBER 99 force field74 are used unless otherwise noted. Periodic boundary conditions are imposed, and the long-range electrostatic and van der Waals interactions are treated with the particle mesh Ewald (PME) method75 and a real-space cutoff at 10 Å, respectively. Bonds involving hydrogen are constrained with the SHAKE algorithm,76 and all simulations (after equilibration) are performed under the constant-NVT (300 K) with a time step of 1 fs, with the temperature maintained with the Langevin thermostat using a collision frequency of 1.0 ps−1. All MD simulations are performed using the AMBER 11 program package,77 with the modifications to implement the IHS method. The Dipole Inversion Model. To examine the efficiency of the weight optimization approaches, we first study the dipole inversion model in explicit water solution as was done in ref 64. To challenge the IHS method, the magnitude of the dipole is substantially enhanced relative to that of ref 64: the diatomic solute has a charge of +1 and −1 on the two atoms, respectively, and they are connected by a bond with a harmonic potential around 1.526 Å with a force constant of 310.0 kcal· mol−1·Å−2; the original model has the partial charges as ±0.25. U0 and U1 are constructed to have identical structures but with

Figure 4. Convergence of weight functions (Ũ i) in the dipole inversion system using the (a) histogram flattening, (b) WHAM, and (c) hybrid approaches. Points in each plot represent where the weight functions are updated; blue, green, red, cyan, purple, and ocher colors represent the weights for λ = 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0, respectively.

To further make use of the samples beyond the adjacent step (i.e., the n-th step), WHAM can be a convenient yet powerful approach. With the use of preceding n samples, the (n + 1)-th weight function can be obtained by iteratively solving the WHAM equations70 n

∑k = 1 N (k){ρi(b)}(k) 1 (n + 1) = − ln Uĩ (l) n β ∑l = 1 N (l) exp[−β(Uĩ − f (l) )] 1 − ln C(n + 1) β

(12)

and N

f

(k)

λ 1 (n + 1) (k) = − ln ∑ exp[−β(Uĩ + Uĩ )] β i=1

(13) 8213



Article

results have a smaller standard deviation by a factor of 6 with less total simulation time. The mean absolute error between the results obtained from IHS and FEP is 0.2 kcal/mol, and that between its symmetric component (with respect to λ = 0.5) is also about 0.2 kcal/mol for both IHS and FEP. QM/MM to MM Perturbation. The ability to quickly obtain the weights is beneficial to the computation of relative free energies of different Hamiltonians, e.g., relative free energy of quantum mechanical (QM) and molecular mechanical (MM) Hamiltonians, as discussed in the Introduction, or pKa calculations29,78,79 (i.e., protonated and deprotonated states). Here we explore this by calculating the relative solvation free energy of QM and MM water molecules in explicit MM water solution. Using the thermodynamic cycle in Figure 1, the relative solvation free energy can be obtained by

opposite charges on the two solute atoms; thus, the free energy as a function of λ is symmetric with respect to λ = 0.5. On the other hand, due to the large difference in the solute dipole upon inversion, a large solvent reorganization is required to switch between the two end states. Thus, this model allows one to evaluate the accuracy of the weight estimate as well as the sampling efficiency over different Uλ minima. Figure 4 compares the time evolution of the weight functions obtained using the (i) histogram flattening, (ii) WHAM, and (iii) hybrid approaches, respectively. In the hybrid approach, the histogram flattening approach is used for the first 150 ps (weight being updated every 1 ps with a cap of 5.0 kcal/mol) and subsequently switched to WHAM, using the last 100 sets of samples to update the weights. The result shows that, while the histogram flattening approach quickly obtains a rough estimate of the weights within only 100 ps, due to the limited sampling between the updates the sample cannot cover the full λ space; thus, the weight fluctuates by a few kcal/mol. On the other hand, WHAM gives very small fluctuations of the weights once they are converged, though due to the poor initial guess (estimated from the first short sampling of 1 ps about only one end state), the optimal value is reached rather slowly for some of the λ states (∼350 ps). In the hybrid approach, the weights are optimized quickly as in the case of (i) histogram flattening, yet the standard deviation becomes small after 200 ps (similar to (ii) WHAM after 450 ps). Once convergence is reached, the amount of sampling per step is extended to obtain a broader distribution over the λ states to further reduce the fluctuation of the weights. For example, in (iii), the sampling time per updating step is increased from 1 to 5 ps after 152 ps and further up to 25 ps after 347 ps. Finally, the estimate of the weights obtained from the last 175 ps of the simulation (8 samples) is used to estimate the free energies and compared to the conventional free energy perturbation (FEP) calculations. In the FEP calculation, 11 λ windows with an increment of 0.1 are used, and each window is sampled for 100 ps (not including the equilibration step of 100 ps each). As shown in Figure 5, the estimates from the two methods are in excellent agreement, and furthermore, the IHS

Although our motivation is to compute such free energy differences involving a high-level ab initio QM/MM potential, for this work, we illustrate the application of IHS by taking two semiempirical QM methods, i.e., AM180 and SCC-DFTB81 (the original second-order formulation82), as the QM part; for the MM model, we explore two popular water models, the TIP3P73 and TIP5P83 models. To simplify the situation, the vdW parameters for the QM atoms are taken to be the same as the target MM model in each perturbation, and the intramolecular geometry of the QM water molecule is fixed to that of the target MM water model as well. Thus, ΔGgas(QM → MM) can be calculated quickly, and these values are used as the initial guesses in the subsequent ΔGaq(QM → MM) calculations using IHS. The time evolutions of the weights for AM1 to TIP3P/ TIP5P transformations (in TIP3P and TIP5P solution, respectively) are shown in Figure 6, and the relative solvation free energies for both AM1 and SCC-DFTB results are

Figure 5. Free energy estimate for λ states in the dipole inversion system. Red and blue bars represent the result from the IHS and conventional free energy perturbation (FEP) calculations, respectively. Free energies are relative to the λ = 0 state, and the error bars are multiplied by 10 to increase visibility. The last 8 samples (175 ps in total) from the IHS result are used to obtain the average, and 100 ps per λ window (excluding 100 ps equilibration) is used for the FEP calculations.

Figure 6. Convergence of weight functions (Ũ i) for QM/MM to MM perturbation of AM1 water into (top) TIP3P and (bottom) TIP5P water models, respectively. Solvent molecules are described by the same MM water model as the solute in each case. Cross points describe when the weights were updated, and blue, green, red, cyan, purple, and ocher lines represent λ = 0.0, 0.2, 0.4, 0.6, 0.8, and 1.0, with U0 and U1 being QM/MM and MM potentials, respectively.

ΔΔGgas → aq (MM → QM)

8214

= ΔGgas → aq (MM) − ΔGgas → aq (QM)

(16)

= ΔGaq (QM → MM) − ΔGgas(QM → MM)

(17)



Article

solution.84 We also see that, compared to AM1, SCC-DFTB gives a higher first peak; especially with the TIP5P water solution, SCC-DFTB shows a clear second solvation shell at around 4.5 Å, in good agreement with the pure TIP5P result. This is likely the reason for the smaller ΔΔG (in absolute value) for SCC-DFTB/TIP5P (3.1 kcal/mol) than SCCDFTB/TIP3P (4.3 kcal/mol) or AM1/MM (>4.5 kcal/mol). Configurational Sampling of Polypeptides. In the third example, IHS is applied to the conformational sampling of a tripeptide Arg-Gly-Asp (RGD) in solution and on a solvated rutile TiO2 (110) surface. This peptide is chosen here because it has been studied in previous work as a sequence that has a high affinity to the TiO2 surface,85−89 although all previous simulations were only up to 4 ns time scale. To enhance the conformational sampling of the RGD peptide, we employ the multistate IHS formulation (eq 15) in which two additional reference Hamiltonians are introduced. In these unphysical reference Hamiltonians, one of the main chain dihedral potentials (ϕ or ψ) for each peptide bond is inverted (by flipping the sign of the force constant), while that for the other (ψ or ϕ) is zeroed out. By design, these unphysical reference Hamiltonians drive the sampling over main chain barriers in the physical system. The force fields for TiO2 are the same as in ref 59, and the same simulation protocol is used to set up the TiO2−RGD complex. In the starting bound configuration, the aspartate side chain adopts the doubly bound mode59 to the Ti atoms on the surface (and this binding mode remains throughout the simulation), since previous studies of amino acid side chain analogues indicated that this is the favorable binding mode of carboxylate and that positively charged residues have weaker binding affinities.59,90 We note that these trends apply because our setup does not consider water ionization at the TiO2 surface and therefore treats the surface as overall charge neutral; this is a simplification91 but sufficient for the purpose of this work. Each IHS simulation is performed for 10 ns; the weights are updated every 1 ps with the histogram flattening method up to 500 ps, and subsequently switched to WHAM using the last 400 steps. We see that the weights converge within the first 100 ps (see Figure S1 in the Supporting Information); thus, the trajectory after 100 ps is considered as the “production” step (while the weights are still updated dynamically). To construct the PMFs, WHAM is applied to the entire production trajectory to determine the relative weight of each snapshot (note that the WHAM protocols for PMF construction and IHS weight function determinations are distinct). As a reference, we perform a 100 ns regular MD sampling in solution and on the surface, respectively. The first 10 ns of these reference trajectories are also used to examine to what degree IHS leads to enhanced sampling. Structures are recorded every 1 ps in the regular simulations, whereas those are extracted every 0.1 ps in IHS. This is because the points from IHS are not equally weighted and more frequent output is desirable to cover the important conformations and to create a smoother PMF surface. The results in solution and on the surface are summarized in Figures 8 and 9, respectively. For the RGD peptide in solution, comparing results for the 10 and 100 ns regular MD simulations (Figure 8b and c) clearly indicates that 10 ns is too short for sampling all important free energy basins. For example, for R and D, the regions of ϕ ∼ 50° are only sampled in the 100 ns simulation; even the ϕ − ψ distributions for the flexible G residue are narrowly sampled at the 10 ns time scale. With 10 ns of IHS (Figure 8a),

summarized in Table 1. Figure 6 shows that in either case the weights converge within 100 ps, with the convergence in the Table 1. Relative Solvation Free Energy ΔΔGgas→aq(QM → MM) (in kcal/mol) with Different Water Models QM method AM1 SCC-DFTB

MM model

ΔGgas

TIP3P TIP5P TIP3P TIP5P

59.2 59.2 247.3 247.3

ΔGaq 54.6 54.4 243.0 244.1

± ± ± ±

ΔΔGgas→aq 0.3 0.0 0.3 0.1

−4.6 −4.8 −4.3 −3.1

± ± ± ±

0.3 0.0 0.3 0.1

TIP5P water case being slightly faster than that in the TIP3P water case. The SCC-DFTB conversion to TIP3P/TIP5P water in solution also shows a similar behavior and thus is not shown. The relative solvation free energies in Table 1 show that, in contrast to the previous ab initio QM/MM → MM results,46 we see negative relative solvation free energies in all cases. This indicates that the semiempirical QM waters studied here are less “soluble” than the MM waters in MM bulk water. The current results do not show a large difference between TIP3P and TIP5P, whereas the ab initio QM/MM results showed a remarkable increase in the relative solvation free energy for TIP5P waters46 (>10 kcal/mol); on the other hand, the latter result might be due, in part, to limited sampling, as indicated from the radial distribution function (RDF) plots in ref 46. To understand the trends of relative solubility of different water models (semiempirical QM < MM < ab initio QM46), we run another 1 ns of IHS simulation using the converged weights. Snapshots are recorded every 0.1 ps, and the RDFs are reconstructed using the reweighting scheme described above (eq 5). The results for all four simulations are shown in Figure 7, which shows that AM1 and SCC-DFTB models have much lower and broader first solvation shells compared to MM (TIP3P and TIP5P) and ab initio QM46 waters. This is consistent with a previous study which reported the tail of the first RDF peak even up to 4 Å for AM1 water in SPC water

Figure 7. Radial distribution function (g(r)) between the “solute” and solvent oxygen atoms from the QM/MM to MM perturbation calculations for the cases of (a) AM1 to TIP3P, (b) AM1 to TIP5P, (c) SCC-DFTB to TIP3P, and (d) SCC-DFTB to TIP5P, respectively. The blue line represents the raw g(r) from the IHS simulation, and the green and red lines represent the reproduced QM/MM and MM results from the IHS result using the reweighting scheme. 8215



Article

Figure 8. Ramchandran plots of R, G, and D (first, second, and third rows, respectively) from the sampling of an RGD peptide in solution: (a) 10 ns of IHS; (b) 10 ns of regular MD; (c) 100 ns of regular MD. Left and right panels describe the sampled points and free energy surface (constructed from the sampled probability), respectively. Colors on the right column range from blue (0 kcal/mol) to red (≥10 kcal/mol). Points on the left panel are plotted every 2 ps.

Figure 9. Ramchandran plots of R, G, and D (first, second, and third rows, respectively) from the sampling of an RGD peptide on a solvated rutile TiO2 (110) surface: (a) 10 ns of IHS; (b) 10 ns of regular MD; (c) 100 ns of regular MD. The left and right panels describe the sampled points and the free energy surface (constructed from the sampled probability), respectively. Colors on the right column range from blue (0 kcal/mol) to red (≥10 kcal/mol). Points on the left panel are plotted every 2 ps. For results for an independent IHS simulation, which has substantially broader sampling for the ϕ, ψ distributions of D, see Figure S2 in the Supporting Information.

around (−60°, −30°) and (50°, 0°) for R and G, respectively. These observations highlight the need of enhanced sampling for even a short peptide at the solid/water interface. Considering the generally encouraging results for IHS for the solution case, the results for the TiO2 surface again suggest that modulating the dihedral potential alone is not sufficient to rapidly sample the conformational changes of a peptide directly bound to an ionic surface; thus, further modulation of peptide/ surface interaction in the reference Hamiltonian in the IHS framework is desirable. This will be explored more systematically in the future. Nevertheless, even without modulating peptide/surface interaction in the reference Hamiltonians, IHS appears to already lead to a substantially better sampling of the peptide conformations than regular MD simulations at the same time scale.

substantially broader distributions are sampled; for example, the ϕ − ψ distribution of G is now very well sampled, and ϕ ∼ 50° regions in R and D plots are both identified as well. On the other hand, the relative free energies (compared to the global minima) of the basins in the IHS results indicate that additional sampling is required to get a quantitative PMF (compared to the 100 ns result). We note that solvent molecules are known to play a role in modulating peptide isomerization dynamics;92,93 thus, more significant enhancement of sampling can likely be achieved if peptide/solvent interactions are also perturbed in the reference Hamiltonians.32,62 For the RGD peptide on the TiO2 surface, the general trends in IHS (Figure 9a) versus regular (Figure 9b and c) MD sampling are similar to the solution case. Regular sampling at the 10 ns time scale misses important low free energy regions, especially for the ϕ − ψ distribution associated with D. The IHS result shown in Figure 9 does not sample the low free energy regions of D with ϕ ∼ −60° either, indicating a bottleneck around ϕ ∼ 0°; an independent IHS trajectory of 10 ns, however, does sample this region (see Figure S2 in the Supporting Information). On the other hand, IHS simulations have located two additional free energy basins that are not sampled at all with the 100 ns regular sampling; the (ϕ, ψ) are

■

CONCLUDING REMARKS Performing adequate sampling is an important aspect of a reliable simulation study of complex systems, and what sampling protocol is effective likely depends on the problem at hand. Motivated by the specific applications of interest (Figure 1) and inspired by the recent work of Gao,61,62 we have developed a novel approach we refer to as integrated 8216



Article

an ab initio QM/MM potential is used. Second, the sampling efficiency of both HREMD and IHS depends on a rapid transition among the different Hamiltonians. In HREMD, this requires a good overlap of distributions among the Hamiltonians and thus potentially a large number of replicas. In IHS, rapid transition among the Hamiltonians is modulated by the weight function, which, as we show using realistic examples (e.g., Figure 4), can be determined on the fly in a robust and automated fashion. In λ-dynamics,24 the coupling parameter that connects the end state Hamiltonians is evolved dynamically during the sampling. When used creatively,94 it can lead to efficient “conformational tunneling” effects. IHS is complementary to λdynamics and related approaches79 in that no fictitious mass needs to be introduced for the λ variable, and transitions among different Hamiltonians are explicitly enhanced by modulating the weight function. In this work, we have only started to explore the applicability of the IHS approach. For practical applications, the efficiency of sampling would clearly depend on the choice of both the “endstate” Hamiltonian(s) (e.g., U0 and U1 in eq 1) and the scheme to construct the intermediate Hamiltonians (i.e., Uλ in eq 1); this requirement is where IHS differs from the original ITS approach,61,62 which does not require specifying any additional Hamiltonian, although one also has to identify the proper subspace in SITS for successful applications.62 For example, for the sampling of a peptide in solution or at the surface, more efficient sampling than reported here can likely be achieved by modulating also the peptide/environment interactions (i.e., a better choice of U1). As for the construction of the intermediate Hamiltonians, we have discussed only a simple linear scheme as used in standard free energy simulations (i.e., eq 1), which may not be most effective at overcoming “hidden barriers” in the “orthogonal space”.25 However, we note inclusion of additional contributions such as the energy gap or other restraints in the IHS framework is possible, leading to sampling protocols complementary to powerful methods5,31 such as orthogonal space random walk25 or umbrella sampling HREMD.95 These developments and a more quantitative comparison to other sampling methods will be explored in future work.

Hamiltonian sampling (IHS), which is complementary to existing techniques for both free energy simulation and enhanced configurational sampling. Using several relatively simple condensed phase examples (Figure 2), we have illustrated the implementation and application of IHS as well as potential developments for the near future. These examples illustrate that, even for seemingly “simple” systems, adequate sampling is essential. For a tri-peptide on the oxide surface, for instance, a regular MD simulation at the 100 ns scale can still miss important configurations. Considering that all previous simulations of this and similar systems that involve oxide surfaces are limited to regular sampling of 5−20 ns, the current observations highlight the need of revisiting some of these problems with more efficient sampling methods. As the name IHS implies, the method constructs an effective Hamiltonian to carry out sampling based on the integrated distributions of a series of Hamiltonians. By judiciously selecting the weights of the different Hamiltonians, one achieves rapid transitions among the energy landscapes that underlie different Hamiltonians and therefore an efficient sampling of important regions of the conformational space. Compared to the original work of (selective) integrated tempering sampling,61,62 IHS takes advantage of distributions of different Hamiltonians rather than distributions at different temperatures; thus, IHS can be used in very flexible ways by choosing different Hamiltonians to drive sampling and provides a complementary technique in applications where very specific sets of degrees of freedom (which can be either conformational or chemical) are to be thoroughly sampled. In this respect, the IHS method shares key features with the enveloping distribution sampling (EDS) approach64,65 pioneered by van Gunsteren and co-workers, who have reported a broad range of successful applications of EDS to biomolecules.63−65,67−69 The main difference between IHS and EDS lies in the way that the distributions of different Hamiltonians are integrated. The parameters in EDS are interdependent and need to be determined by a single long reference simulation. By contrast, the weights in IHS are independent (other than a normalization condition) and can be determined on the f ly, with a well-defined criterion (e.g., eq 6), by integrating multiple trajectories via WHAM or MBAR protocols. As a result, introduction of multiple Hamiltonians in IHS for additional flexibility and efficiency is straightforward. Nevertheless, the key motivations for IHS and EDS are very similar. The IHS approach can also be compared to several other related sampling methods, especially HREMD and λ-dynamics. In HREMD,14,15 different Hamiltonians are used to drive the sampling of different regions of the conformational space; thus, like IHS, the efficiency depends on the choice of the Hamiltonians. The key difference is that HREMD involves running many parallel simulations, each governed by a different Hamiltonian, while IHS runs a single simulation with an effective Hamiltonian determined by the integrated distributions of the end-state and intermediate Hamiltonians. This difference has two consequences. First, the computational cost for HREMD clearly increases linearly as the number of replicas, while, in IHS, if the intermediate Hamiltonians are a simple algebraic combination of the end states (e.g., eq 1), the computational cost varies little with respect to the number of intermediate Hamiltonians. Although this saving is not compelling when large computational resources are available for parallel computations, it can still be advantageous when the potential functions are computationally expensive, such as when

■

ASSOCIATED CONTENT

* Supporting Information S

Convergence behavior of the weights in the RGD peptide simulations is illustrated using the solution case. Results for an independent IHS simulation for the RGD peptide at the solvated TiO2 surface are also included. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Authors

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

Institute for Molecular Science, Myodaiji, Okazaki, Aichi, 4448585, Japan. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is dedicated to Professor Jim Skinner, a much appreciated and admired colleague, in honor of his numerous 8217



Article

contributions to the field of statistical mechanics and his leadership as the director of the Theoretical Chemistry Institute at the University of Wisconsin, Madison. T.M. appreciates the support of the JSPS Postdoctoral Fellowship for Research Abroad. The method development is partially funded by a National Science Foundation (NSF) grant CHE-1300209 to Q.C., and the application to model adsorption on rutile TiO2 surfaces is supported by NSF-ECS-1152604 to J.A.P. and R.J.H. Computational resources from the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF grant number OCI-1053575, are greatly appreciated; computations are also supported in part by NSF through a major instrumentation grant (CHE-0840494) to the Chemistry department.

■

(19) Maragliano, L.; vanden-Eijnden, E. A temperature accelerated method for sampling free energy and determining reaction pathways in rare events simulations. Chem. Phys. Lett. 2006, 426, 168−175. (20) Huang, X. H.; Bowman, G. R.; Bacallado, S.; Pande, V. S. Rapid equilibrium sampling initiated from nonequilibrium data. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 19765−19769. (21) Zheng, W. W.; Rohrdanz, M. A.; Clementi, C. Rapid Exploration of Configuration Space with Diffusion-Map-Directed Molecular Dynamics. J. Phys. Chem. B 2013, 117, 12769−12776. (22) Kim, J.; Straub, J. E.; Keyes, T. Statistical-temperature monte carlo and molecular dynamics algorithms. Phys. Rev. Lett. 2006, 97, 050601. (23) Tian, P.; Jónsson, S. A.; Ferkinghoff-Borg, J.; Krivov, S. V.; Lindorff-Larsen, K.; Irbäck, A.; Boomsma, W. Robust estimation of diffusion-optimized ensembles for enhanced sampling. J. Chem. Theory Comput. 2014, 10, 543−553. (24) Kong, X.; Brooks, C. L. λ-dynamics: A new approach to free energy calculations. J. Chem. Phys. 1996, 105, 2414. (25) Zheng, L.; Chen, M.; Yang, W. Random walk in orthogonal space to achieve efficient free-energy simulation of complex systems. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 20227−20232. (26) Woo, H.-J.; Roux, B. Calculation of absolute protein-ligand binding free energy from computer simulations. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6825−6830. (27) Henin, J.; Fiorin, G.; Chipot, C.; Klein, M. L. Exploring Multidimensional Free Energy Landscapes Using Time-Dependent Biases on Collective Variables. J. Chem. Theory Comput. 2010, 6, 35− 47. (28) Pohorille, A.; Jarzynski, C.; Chipot, C. Good Practices in FreeEnergy Calculations. J. Phys. Chem. B 2010, 114, 10235−10253. (29) Itoh, S. G.; Damjanović, A.; Brooks, B. R. pH replica-exchange method based on discrete protonation states. Proteins 2011, 79, 3420− 3436. (30) Hu, Y.; Hong, W.; Shi, Y.; Liu, H. Temperature-Accelerated Sampling and Amplified Collective Motion with Adiabatic Reweighting to Obtain Canonical Distributions and Ensemble Averages. J. Chem. Theory Comput. 2012, 8, 3777−3792. (31) Chipot, C. Frontiers in free-energy calculations of biological systems. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4, 71−89. (32) Liu, P.; Kim, B.; Friesner, R. A.; Berne, B. J. Replica exchange with solute tempering: a method for sampling biological systems in explicit water. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 13749−13754. (33) Bolhuis, P. G. Rare events via multiple reaction channels sampled by path replica exchange. J. Chem. Phys. 2008, 129, 114108. (34) Jiang, W.; Roux, B. Free Energy Perturbation Hamiltonian Replica-Exchange Molecular Dynamics (FEP/H-REMD) for Absolute Ligand Binding Free Energy Calculations. J. Chem. Theory Comput. 2010, 6, 2559−2565. (35) Kokubo, H.; Tanaka, T.; Okamoto, Y. Two-Dimensional Replica-Exchange Method for Predicting Protein-Ligand Binding Structures. J. Comput. Chem. 2013, 34, 2601−2614. (36) Mitsutake, A.; Okamoto, Y. Replica-exchange simulated tempering method for simulations of frustrated systems. Chem. Phys. Lett. 2000, 332, 131−138. (37) Zhang, Y. Pseudobond ab initio QM/MM approach and its applications to enzyme reactions. Theor. Chem. Acc. 2006, 116, 43−50. (38) Hu, H.; Yang, W. T. Free Energies of Chemical Reactions in Solution and in Enzymes with Ab Initio Quantum Mechanics/ Molecular Mechanics Methods. Annu. Rev. Phys. Chem. 2008, 59, 573− 601. (39) Kamerlin, S. C. L.; Haranczyk, M.; Warshel, A. Progress in Ab Initio QM/MM Free-Energy Simulations of Electrostatic Energies in Proteins: Accelerated QM/MM Studies of pK(a), Redox Reactions and Solvation Free Energies. J. Phys. Chem. B 2009, 113, 1253−1272. (40) Senn, H. M.; Thiel, W. QM/MM methods for biomolecular systems. Angew. Chem., Int. Ed. 2009, 48, 1198−1229. (41) Laio, A.; VandeVondele, J.; Rothlisberger, U. A Hamiltonian electrostatic coupling scheme for hybrid Car-Parrinello molecular dynamics simulations. J. Chem. Phys. 2002, 116, 6941−6947.

REFERENCES

(1) Kollman, P. Free energy calculations: Applications to chemical and biochemical phenomena. Chem. Rev. 1993, 93, 2395−2417. (2) Karplus, M.; Kuriyan, J. Molecular dynamics and protein function. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6679−6685. (3) van Gunsteren, W. F.; et al. Biomolecular modeling: Goals, problems, perspectives. Angew. Chem., Int. Ed. 2006, 45, 4064−4092. (4) Adcock, S. A.; McCammon, J. A. Molecular dynamics: Survey of methods for simulating the activity of proteins. Chem. Rev. 2006, 106, 1589−1615. (5) Zuckerman, D. M. Equilibrium Sampling in Biomolecular Simulations. Annu. Rev. Biophys. 2011, 40, 41−62. (6) Dror, R. O.; Dirks, R. M.; Grossman, J. P.; Xu, H.; Shaw, D. E. Biomolecular simulation: a computational microscope for molecular biology. Annu. Rev. Biophys. 2012, 41, 429−452. (7) Huang, X.; Hagen, M.; Kim, B.; Friesner, R. A.; Zhou, R.; Berne, B. J. Replica Exchange with Solute Tempering: Efficiency in Large Scale Systems. J. Phys. Chem. B 2007, 111, 5405−5410. (8) Hamelberg, D.; Mongan, J.; McCammon, J. Accelerated molecular dynamics: A promising and efficient simulation method for biomolecules. J. Chem. Phys. 2004, 120, 11919. (9) Park, S.; Schulten, K. Calculating potentials of mean force from steered molecular dynamics simulations. J. Chem. Phys. 2004, 120, 5946. (10) Marsili, S.; Barducci, A.; Chelli, R.; Procacci, P.; Schettino, V. Self-healing Umbrella Sampling: A Non-equilibrium Approach for Quantitative Free Energy Calculations. J. Phys. Chem. B 2006, 110, 14011−14013. (11) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562−12566. (12) Zhang, C.; Ma, J. Enhanced sampling and applications in protein folding in explicit solvent. J. Chem. Phys. 2010, 132, 244101. (13) Zheng, L.; Yang, W. Essential energy space random walks to accelerate molecular dynamics simulations: Convergence improvements via an adaptive-length self-healing strategy. J. Chem. Phys. 2008, 129, 014105. (14) Sugita, Y.; Kitao, A.; Okamoto, Y. Multidimensional replicaexchange method for free-energy calculations. J. Chem. Phys. 2000, 113, 6042. (15) Fukunishi, H.; Watanabe, O.; Takada, S. On the Hamiltonian replica exchange method for efficient sampling of biomolecular systems: Application to protein structure prediction. J. Chem. Phys. 2002, 116, 9058−9067. (16) Lyubartsev, A. P.; Martsinovski, A. A.; Shevkunov, S. V.; Vorontsov-Velyaminov, P. N. New approach to Monte Carlo calculation of the free energy: Method of expanded ensembles. J. Chem. Phys. 1992, 96, 1776. (17) Marinari, E.; Parisi, G. Simulated Tempering: A New Monte Carlo Scheme. Europhys. Lett. 2007, 19, 451−458. (18) Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 1999, 314, 141−151. 8218



Article

(42) Field, M. J.; Bash, P. A.; Karplus, M. A Combined QuantumMechanical and Molecular Mechanical Potential for MolecularDynamics Simulations. J. Comput. Chem. 1990, 11, 700−733. (43) Gao, J. In Reviews in Computational Chemistry VII; Lipkowitz, K. B., Boyd, D. B., Eds.; VCH: New York, 1995; p 119. (44) Monard, G.; Merz, K. M., Jr. Combined quantum mechanical/ molecular mechanical methodologies applied to biomolecular systems. Acc. Chem. Res. 1999, 32, 904−911. (45) Riccardi, D.; Schaefer, P.; Yang, Y.; Yu, H.; Ghosh, N.; PratResina, X.; Konig, P.; Li, G.; Xu, D.; Guo, H.; Elstner, M.; Cui, Q. Development of effective quantum mechanical/molecular mechanical (QM/MM) methods for complex biological processes (Feature Article). J. Phys. Chem. B 2006, 110, 6458−6469. (46) Shaw, K. E.; Woods, C. J.; Mulholland, A. J. Compatibility of Quantum Chemical Methods and Empirical (MM) Water Models in Quantum Mechanics/Molecular Mechanics Liquid Water Simulations. J. Phys. Chem. Lett. 2010, 1, 219−223. (47) König, G.; Hudson, P. S.; Boresch, S.; Woodcock, H. L. Multiscale free energy simulations: An efficient method for connecting classical MD simulations to QM or QM/MM free energies using NonBoltzmann Bennett reweighting schemes. J. Chem. Theory Comput. 2014, DOI: 10.1021/ct401118k. (48) Schravendijk, P.; Ghiringhelli, L. M.; Site, L. D.; van der Vegt, N. F. Interaction of hydrated amino acids with metal surfaces: A multiscale modeling description. J. Phys. Chem. C 2007, 111, 2631− 2642. (49) Ghiringhelli, L. M.; Hess, B.; van der Vegt, N. F. A.; Delle Site, L. Competing adsorption between hydrated peptides and water onto metal surfaces: from electronic to conformational properties. J. Am. Chem. Soc. 2008, 130, 13460−13464. (50) Di Felice, R.; Corni, S. Simulation of Peptide-Surface Recognition. J. Phys. Chem. Lett. 2011, 2, 1510−1519. (51) Lower, B. H.; Lins, R. D.; Oestreicher, Z.; Straatsma, T. P.; Hochella, M. F.; Shi, L.; Lower, S. K. In vitro evolution of a peptide with a hematite binding motif that may constitute a natural metaloxide binding archetype. Environ. Sci. Technol. 2008, 42, 3821−3827. (52) Goede, K.; Busch, P.; Grundmann, M. Binding Specificity of a Peptide on Semiconductor Surfaces. Nano Lett. 2004, 4, 2115−2120. (53) Hayashi, T.; Sano, K.-I.; Shiba, K.; Kumashiro, Y.; Iwahori, K.; Yamashita, I.; Hara, M. Mechanism underlying specificity of proteins targeting inorganic materials. Nano Lett. 2006, 6, 515−519. (54) Imamura, K.; Kawasaki, Y.; Nagayasu, T.; Sakiyama, T.; Nakanishi, K. Adsorption characteristics of oligopeptides composed of acidic and basic amino acids on titanium surface. J. Biosci. Bioeng. 2007, 103, 7−12. (55) Shiba, K. Exploitation of peptide motif sequences and their use in nanobiotechnology. Curr. Opin. Biotechnol. 2010, 21, 412−425. (56) Sano, K.-I.; Shiba, K. A Hexapeptide Motif that Electrostatically Binds to the Surface of Titanium. J. Am. Chem. Soc. 2003, 125, 14234− 14235. (57) Naik, R. R.; Brott, L. L.; Clarson, S. J.; Stone, M. O. Silicaprecipitating peptides isolated from a combinatorial phage display peptide library. J. Nanosci. Nanotechnol. 2002, 2, 95−100. (58) Patwardhan, S. V.; Emami, F. S.; Berry, R. J.; Jones, S. E.; Naik, R. R.; Deschaume, O.; Heinz, H.; Perry, C. C. Chemistry of Aqueous Silica Nanoparticle Surfaces and the Mechanism of Selective Peptide Adsorption. J. Am. Chem. Soc. 2012, 134, 6244−6256. (59) Mori, T.; Hamers, R. J.; Pedersen, J. A.; Cui, Q. An Explicit Consideration of Desolvation is Critical to Binding Free Energy Calculations of Charged Molecules at Ionic Surfaces. J. Chem. Theory Comput. 2013, 9, 5059−5069. (60) Yang, L. J.; Shao, Q.; Gao, Y. Q. Comparison between integrated and parallel tempering methods in enhanced sampling simulations. J. Chem. Phys. 2009, 130, 124111. (61) Gao, Y. Q. An integrate-over-temperature approach for enhanced sampling. J. Chem. Phys. 2008, 128, 064105. (62) Yang, L.; Gao, Y. Q. A selective integrated tempering method. J. Chem. Phys. 2009, 131, 214109.

(63) Christ, C. D.; van Gunsteren, W. F. Multiple free energies from a single simulation: Extending enveloping distribution sampling to nonoverlapping phase-space distributions. J. Chem. Phys. 2008, 128, 174112. (64) Christ, C. D.; van Gunsteren, W. F. Enveloping distribution sampling: A method to calculate free energy differences from a single simulation. J. Chem. Phys. 2007, 126, 184110. (65) Christ, C. D.; van Gunsteren, W. F. Simple, Efficient, and Reliable Computation of Multiple Free Energy Differences from a Single Simulation: A Reference Hamiltonian Parameter Update Scheme for Enveloping Distribution Sampling (EDS). J. Chem. Theory Comput. 2009, 5, 276−286. (66) Oostenbrink, C.; van Gunsteren, W. F. Free energies of ligand binding for structurally diverse compounds. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6750−6754. (67) Christ, C. D.; van Gunsteren, W. F. Comparison of three enveloping distribution sampling Hamiltonians for the estimation of multiple free energy differences from a single simulation. J. Comput. Chem. 2009, 30, 1664−1679. (68) Lin, Z. X.; van Gunsteren, W. F. Enhanced conformational sampling using enveloping distribution sampling. J. Chem. Phys. 2013, 139, 144105. (69) van Gunsteren, W. F. Combination of Enveloping Distribution Sampling (EDS) of a Soft-Core Reference-State Hamiltonian with One-Step Perturbation to Predict the Effect of Side Chain Substitution on the Relative Stability of Right- and Left-Helical Folds of betaPeptides. J. Chem. Theory Comput. 2013, 9, 126−134. (70) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. THE weighted histogram analysis method for freeenergy calculations on biomolecules. I. The method. J. Comput. Chem. 1992, 13, 1011−1021. (71) Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. J. Chem. Phys. 2008, 129, 124105. (72) Wang, F.; Landau, D. Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States. Phys. Rev. Lett. 2001, 86, 2050−2053. (73) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (74) Cornell, W.; Cieplak, P.; Bayly, C.; Gould, I.; Merz, K.; Ferguson, D.; Spellmeyer, D.; Fox, T.; Caldwell, J.; Kollman, P. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules (vol 117, pg 5179, 1995). J. Am. Chem. Soc. 1996, 118, 2309−2309. (75) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577−8593. (76) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327−341. (77) Case, D. A.; et al. AMBER 12; 2012. (78) Riccardi, D.; Schaefer, P.; Cui, Q. pKa calculations in solution and proteins with QM/MM free energy perturbation simulations. J. Phys. Chem. B 2005, 109, 17715−17733. (79) Donnini, S.; Tegeler, F.; Groenhof, G.; Grubmüller, H. Constant pH Molecular Dynamics in Explicit Solvent with lambdaDynamics. J. Chem. Theory Comput. 2011, 7, 1962−1978. (80) Thiel, W. Perspectives on semiempirical molecular orbital theory. Adv. Chem. Phys. 1996, 93, 703−757. (81) Gaus, M.; Cui, Q.; Elstner, M. Density Functional Tight Binding (DFTB): Application to organic and biological molecules. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4, 49−61. (82) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-consistent-charge densityfunctional tight-binding method for simulations of complex materials properties. Phys. Rev. B 1998, 58, 7260−7268. 8219



Article

(83) Mahoney, M. W.; Jorgensen, W. L. A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions. J. Chem. Phys. 2000, 112, 8910. (84) Geerke, D. P.; Thiel, S.; Thiel, W.; van Gunsteren, W. F. QMMM interactions in simulations of liquid water using combined semiempirical/classical Hamiltonians. Phys. Chem. Chem. Phys. 2007, 10, 297. (85) Rammelt, S.; Illert, T.; Bierbaum, S.; Scharnweber, D.; Zwipp, H.; Schneiders, W. Coating of titanium implants with collagen, RGD peptide and chondroitin sulfate. Biomaterials 2006, 27, 5561−5571. (86) Zhang, H.-p.; Lu, X.; Leng, Y.; Watari, F.; Weng, J.; Feng, B.; Qu, S. Effects of aqueous environment and surface defects on Arg-GlyAsp peptide adsorption on titanium oxide surfaces investigated by molecular dynamics simulation. J. Biomed. Mater. Res. 2010, 96A, 466− 476. (87) Chen, M.; Wu, C.; Song, D.; Li, K. RGD tripeptide onto perfect and grooved rutile surfaces in aqueous solution: adsorption behaviors and dynamics. Phys. Chem. Chem. Phys. 2009, 12, 406−415. (88) Schneider, J.; Ciacchi, L. C. A Classical Potential to Model the Adsorption of Biological Molecules on Oxidized Titanium Surfaces. J. Chem. Theory Comput. 2011, 7, 473−484. (89) Song, D.-P.; Chen, M.-J.; Liang, Y.-C.; Bai, Q.-S.; Chen, J.-X.; Zheng, X.-F. Adsorption of tripeptide RGD on rutile TiO2 nanotopography surface in aqueous solution. Acta Biomater. 2010, 6, 684−694. (90) Monti, S.; Walsh, T. R. Free Energy Calculations of the Adsorption of Amino Acid Analogues at the Aqueous Titania Interface. J. Phys. Chem. C 2010, 114, 22197−22206. (91) Cheng, J.; Sprik, M. Acidity of the Aqueous Rutile TiO2 (110) Surface from Density Functional Theory Based Molecular Dynamics. J. Chem. Theory Comput. 2010, 6, 880−889. (92) Bolhuis, P.; Dellago, C.; Chandler, D. Reaction coordinates of biomolecular isomerization. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 5877−5882. (93) Ma, A.; Dinner, A. R. Automatic method for identifying reaction coordinates in complex systems. J. Phys. Chem. B 2005, 109, 6769− 6779. (94) Li, H. Z.; Fajer, M.; Yang, W. Simulated scaling method for localized enhanced sampling and simultaneous “alchemical” free energy simulations: A general method for molecular mechanical, quantum mechanical, and quantum mechanical/molecular mechanical simulations. J. Chem. Phys. 2007, 126, 024106. (95) Murata, K.; Sugita, Y.; Okamoto, Y. Free energy calculations for DNA base stacking by replica-exchange umbrella sampling. Chem. Phys. Lett. 2004, 385, 1−7.

8220


Integrated Hamiltonian Sampling: A Simple and Versatile Method for

Recommend Documents