Predictive Sampling of Rare Conformational Events in Aqueous

Article pubs.acs.org/JCTC

Predictive Sampling of Rare Conformational Events in Aqueous Solution: Designing a Generalized Orthogonal Space Tempering Method Chao Lu,†,§,∥ Xubin Li,†,§ Dongsheng Wu,‡ Lianqing Zheng,*,‡ and Wei Yang*,†,‡ †

Department of Chemistry and Biochemistry, Florida State University, Tallahassee, Florida 32306, United States Institute of Molecular Biophysics, Florida State University, Tallahassee, Florida 32306, United States

‡

ABSTRACT: In aqueous solution, solute conformational transitions are governed by intimate interplays of the fluctuations of solute−solute, solute−water, and water−water interactions. To promote molecular fluctuations to enhance sampling of essential conformational changes, a common strategy is to construct an expanded Hamiltonian through a series of Hamiltonian perturbations and thereby broaden the distribution of certain interactions of focus. Due to a lack of active sampling of configuration response to Hamiltonian transitions, it is challenging for common expanded Hamiltonian methods to robustly explore solvent mediated rare conformational events. The orthogonal space sampling (OSS) scheme, as exemplified by the orthogonal space random walk and orthogonal space tempering methods, provides a general framework for synchronous acceleration of slow configuration responses. To more effectively sample conformational transitions in aqueous solution, in this work, we devised a generalized orthogonal space tempering (gOST) algorithm. Specifically, in the Hamiltonian perturbation part, a solvent-accessible-surfacearea-dependent term is introduced to implicitly perturb near-solute water−water fluctuations; more importantly in the orthogonal space response part, the generalized force order parameter is generalized as a two-dimension order parameter set, in which essential solute−solvent and solute−solute components are separately treated. The gOST algorithm is evaluated through a molecular dynamics simulation study on the explicitly solvated deca-alanine (Ala10) peptide. On the basis of a fully automated sampling protocol, the gOST simulation enabled repetitive folding and unfolding of the solvated peptide within a single continuous trajectory and allowed for detailed constructions of Ala10 folding/unfolding free energy surfaces. The gOST result reveals that solvent cooperative fluctuations play a pivotal role in Ala10 folding/unfolding transitions. In addition, our assessment analysis suggests that because essential conformational events are mainly driven by the compensating fluctuations of essential solute−solvent and solute−solute interactions, commonly employed “predictive” sampling methods are unlikely to be effective on this seemingly “simple” system. The gOST development presented in this paper illustrates how to employ the OSS scheme for physics-based sampling method designs.

■

(TREM)11,12 and Hamiltonian replica exchange methods (HREM)13−16] and expanded ensemble methods [simulated tempering (ST) method17−19 and simulated scaling (SS) method20] etc., thereby to broaden the distribution of certain interactions of focus. Taking the potential scaling scheme, such as represented by HREM and SS, as the example, the general expanded Hamiltonian can be expressed as

INTRODUCTION Governed by the intimate interplays of solute−solute, solute− water, and water−water interactions, solute molecules in aqueous solution constantly fluctuate. Critical fluctuations, when reaching adequate amplitudes, lead to significant conformational transitions. Understanding conformational changes in aqueous solution has been a key topic in the field of biomolecular simulation. For the fact that many biomolecular processes are beyond time scales that canonical ensemble molecular dynamics (MD) can routinely reach, specialized sampling techniques are of great importance. The past decades have observed the flourishing of enhancedsampling method developments.1−7 An important aspect of these efforts is to develop “predictive” sampling methods that do not rely on system-specific knowledge. A general strategy is to construct an expanded Hamiltonian through a series of perturbed Hamiltonians, as is exemplified by replica exchange8−10 methods [temperature replica exchange method © 2015 American Chemical Society

Hλ = λUss(XN ) + SE(λ)Use(XN ) + Ue(XN ) N

+

∑ Pi 2 /2mi + pλ 2 /2mλ i=1

(1)

where Uss(XN) and Use(XN) respectively represent selected solute−solute and solute−solvent interactions that are subject Received: October 7, 2015 Published: November 12, 2015 41

DOI: 10.1021/acs.jctc.5b00953 J. Chem. Theory Comput. 2016, 12, 41−52

Article

Journal of Chemical Theory and Computation

eq 2 is specifically defined as Uss(XN) + [dSE(λ)/dλ]Use(XN); it has no direct solvent−solvent contribution but two solute related components, between which the ratio is strictly SE(λ) dependent. As discussed in our early study,32 along different microscopic events, the relative fluctuations between Uss(XN) and Use(XN) may follow drastically distinct directions. The actual relative fluctuation directions may not be in good accord with their ratio [1:dSE(λ)/dλ] in the response order parameter function Fλ. Therefore, through the original OST treatment, efficient activation of slow essential configuration responses cannot be guaranteed. In this work, to more robustly accelerate configuration response sampling in aqueous environments, we devised a generalized orthogonal space tempering (gOST) algorithm. Specifically, in the Hamiltonian perturbation part, besides a common potential scaling treatment, a solventaccessible-surface-area-dependent term is introduced to implicitly perturb near-solute water−water fluctuations; more importantly in the orthogonal space response part, the response order parameter is generalized as a two-dimension order parameter set, where the SE(λ)-scaled solute−solvent interaction is separately treated. Through these generalization treatments, first, water−water response fluctuations can be actively accelerated and, second, irrespective of function form chosen for SE(λ), response fluctuations of essential solute− solute and solute−solvent interactions can be robustly enlarged. The gOST algorithm was evaluated through a molecular dynamics simulation study on the explicitly solvated decaalanine (Ala10) peptide. Based on a fully automated sampling protocol, the gOST simulation enabled repetitive folding and unfolding of the solvated peptide and allowed for detailed constructions of Ala10 folding/unfolding free energy surfaces. The gOST result reveals that solvent cooperative fluctuations play a pivotal role in Ala10 folding/unfolding transitions. In addition, our assessment analysis suggests that because essential rare events are mainly driven by compensating fluctuations of essential solute−solvent and solute−solute interactions, common sampling methods are unlikely to be effective on this seemingly simple system. The gOST development presented in this paper illustrates the advantage of the OSS scheme for physics-based predictive sampling method designs.

to the scaling perturbation treatment and Ue(XN) represents the remaining “environment” interactions. To be general, while Uss(XN) is straightly scaled by a factor λ, Use(XN) is scaled by a prechosen λ-dependent function SE(λ). The λ-space dynamics can be realized through either hybrid Monte Carlo (HMC)21 based methods or the extended dynamics approach,22−24 the latter of which requires the scaling factor λ to be treated as a one-dimension particle with its kinetic energy of (Pλ2/2mλ). To enable random and broad visits of prechosen λ states, one usually needs to introduce a biasing potential function f m(λ) (for the HMC scheme, correspondingly a biasing weight function) or utilize the MC-based replica exchange technique.8−10 In comparison with “global” sampling methods like TREM and ST, such potential-scaling-based “focusing” sampling approaches allow fluctuation promotions to be selectively concentrated on essential interactions and generate much less diffusion-sampling overheads. In the past decades, the potential-scaling-based expanded Hamiltonian strategy has attracted tremendous interest. As is increasingly known, applying these approaches to explore rare conformational events in aqueous solution is not readily trivial, even on systems with moderate complexities. Indeed, although varying λ may lead to a shifted distribution, it can take very long time for the system to adequately evolve so as to realize the distribution change; i.e. slow configuration responses (to λ transitions), particularly ones involving rare activated transitions, are likely to be the major sampling bottleneck. To tackle this configuration “lagging” issue, which in the context of free energy simulation is also called the “hidden free energy barrier” problem,25 we introduced the orthogonal space sampling (OSS) scheme,25−27 as exemplified by the orthogonal space random walk (OSRW)25,26 and orthogonal space tempering (OST)27 methods. To enable OST sampling in the context of the potential-scaling Hamiltonian simulation, the Hamiltonian in eq 1 needs to be further modified as the following: Hm = Hλ + fm (λ) + αgm(λ , Fλ)

(2)

where to achieve λ space random walks, the target function of f m(λ) is set as −Go(λ), the negative of the λ-dependent free energy function under the Hλ ensemble, and to synchronously accelerate configuration responses, the target function of gm(λ,Fλ) is set as −Gf(λ,Fλ), the negative of the free energy function along (λ,Fλ) under the Hλ − Go(λ) ensemble. In OST, the generalized force Fλ, specifically defined as ∂Hλ/∂λ, plays the role of the response order parameter. Here, the parameter α controls the magnitude of response fluctuation enlargement; when α is equal to one, the OST method turns to the OSRW method. It is worth noting that α can be expressed as (TES − To)/TES, where To is the system reservoir temperature and TES is called the orthogonal space sampling temperature27 because at each λ state of the Hm ensemble, the Fλ distribution is proportional to e−Gf(λ,Fλ)/(kTES). As was previously noted,25 the original OSS treatment can only robustly accelerate strongly coupled response fluctuations. In the context of the expanded Hamiltonian sampling, its performance is highly sensitive to the specific form of the scaling function SE(λ), which unfortunately needs to be empirically chosen; for instance, with the general condition of “SE(1) = 1,” commonly employed functional forms include λ, (λ + 1)/2,15 λ1/2,28−31 etc. As one can tell from eq 1, SE(λ) governs how different interactions, including unscaled solvent− solvent interactions, dynamically interplay at different perturbed states. More to note, the response order parameter Fλ in

■

THEORETICAL DESIGN In this section, we first summarize the original orthogonal space tempering (OST) method, then introduce the new generalized orthogonal space tempering (gOST) algorithm, and finally describe how to construct free energy surfaces based on samples collected from gOST simulations. Brief Introduction to the Orthogonal Space Tempering (OST) Method. The OST algorithm is an adaptive sampling method. To recursively generate the target biasing energy functions, an extended-dynamics based “doubleintegration” recursion strategy was developed.27 To do so, the target OST Hamiltonian in eq 2 needs to be modified as follows: Hm = Hλ +

pϕ 2 2mϕ

+

+ αgm(λ , ϕ)

⎞2 1 ⎛ ∂Hλ − ϕ⎟ + fm (λ) Kϕ⎜ ⎠ 2 ⎝ ∂λ (3)

in which, by introducing an additional fictitious particle ϕ, we extend the expanded Hamiltonian Hλ in eq 1 to Hϕ = Hλ + pϕ2/ 2mϕ + (1/2)Kϕ(∂Hλ/∂λ − ϕ)2; here through the restraint term, 42


Article

Journal of Chemical Theory and Computation (1/2)Kϕ(∂Hλ/∂λ − ϕ)2, ϕ can play a role as the dynamic restraining reference (DRR) of ∂Hλ/∂λ. The same as the scaling factor λ, the DDR variable ϕ needs to be treated as a onedimension particle with a mass of mϕ. In OST, the extended particles, λ and ϕ, are propagated based on independent Langevin dynamics equations, the temperatures in which should be set as the same as the system reservoir temperature To. To realize the “double-integration” recursion, for each λ′ state, gm(λ′,ϕ) can be adaptively evaluated via an adaptive biasing force (ABF)33 like procedure, specifically based on the following equation:

dλ)Usasa(XN)] to a two-dimension response order parameter set inter intra (Fintra denotes Uss(XN) + (dSASA(λ)/dλ) λ , Fλ ), where Fλ inter Usasa(XN), and Fλ denotes (dSE(λ)/dλ)Use(XN). As discussed in our previous study,32 following the progress of microscopic events, solute−solute and solute−solvent interactions may have distinct interplays and often fluctuate in a compensating manner; therefore separately treating them through a twodimension order parameter set can lead to more robust activation of their correlated fluctuations. Generalized from eq 2, the gOST Hamiltonian should be expressed as

∂Hϕ

in which Hλ is the same as in eq 7 and the target function of inter intra inter gm(λ,Fintra Fλ ), the negative of the λ ,Fλ ) is set as −Gf(λ,Fλ intra inter free energy function along (λ,Fλ ,Fλ ) under the ensemble of Hλ − G0(λ). To realize the double-integration recursion, the gOST Hamiltonian needs to be further modified as

gm(λ′, ϕ) = −

∫ϕ

ϕ

dϕ + Cg (λ′)

∂ϕ

min(λ′)

Hm = Hλ + fm (λ) + αgm(λ , Fλintra , Fλinter)

m

λ′

(4)

where the integration is carried out within the dynamically evolving boundary [ϕmin(λ′),ϕmax(λ′)]; and the calibration constant Cgm(λ′) is calculated as follows: Cg (λ′) = − max[gm(λ′, ϕ)]

Hm = Hλ +

m

− kTo ln

∫ϕ

ϕmax (λ′)

min(λ′)

⎛ g (λ′, ϕ) − max[g (λ′, ϕ)] ⎞ m exp⎜⎜ m ⎟⎟ dϕ kTo ⎝ ⎠

+ (5)

fm (λ) = −

∫ λmin

ϕ (λ′)

( exp(

∫ϕ max(λ′) ϕ exp min

ϕ (λ′)

∫ϕ max(λ′) min

gm(λ′, ϕ) kTo

) dϕ dλ′ ) dϕ

+

2m ϕ1

(6)

+

pϕ 2 1

2m ϕ1

2m ϕ2

1 K ϕ (Fλintra − ϕ1)2 2 1

(9)

+

pϕ 2 2

2m ϕ2

+

1 K ϕ (Fλintra − ϕ1)2 2 1

1 K ϕ (Fλinter − ϕ2)2 2 2

N i=1

+

Different from ϕ in the original OST recursion, which serves as a dynamic restraining reference, here, ϕ1 and ϕ2 serve as adaptive dynamic reporters (ADRs). Specifically, during each recursion, we adaptively generate gm(λ,ϕ1,ϕ2), whereas during each energy calculation and dynamic propagation, we use Fintra λ and Finter instead of ϕ1 and ϕ2 as the respective variables. As for λ the “double-integration” recursion, for each λ′ state, we can obtain gm(λ′,ϕ1,ϕ2) via the following path integration:

Hλ = λUss(XN ) + SE(λ)Use(XN ) + [SASA(λ) − 1]

∑ Pi 2 /2mi + pλ 2 /2mλ

2

1 K ϕ (Fλinter − ϕ2)2 + fm (λ) 2 2

Hϕ1, ϕ2 = Hλ +

It is noted that for OST free energy calculations, an extra correction operation is required to obtain accurate Go(λ); for predictive sampling, the above equation should be adequate. Generalized Orthogonal Space Tempering (gOST) Method. As discussed earlier, in gOST, two major changes are made from the original OST method. First, a solventaccessible-surface-area-dependent term, [SASA(λ)−1] Usasa(XN), is introduced to construct a new expanded Hamiltonian:

Usasa(XN ) + Ue(XN ) +

pϕ 2

in which, instead of ϕ, we introduced two extended particles, ϕ1 and ϕ2. Then, the Hλ Hamiltonian in eq 7 is generalized to a new Hamiltonian:

gm(λ′, ϕ) kTo

1

+ αgm(λ , Fλintra , Fλinter)

Then the two-dimension function gm(λ,ϕ) can be generated via the assembling of calculated one-dimension gm(λ′,ϕ) functions. Furthermore, we can adaptively acquire f m(λ) through another integration calculation: λ

pϕ 2

(8)

gm(λ′, ϕ1 , ϕ2) = − (7)

1,min(λ′), ϕ2,min(λ′)]

+ Cg (λ′)

where Usasa(XN) represents the solvent accessible surface area (SASA) energy that is commonly used to model nonpolar solvation contributions34 and SASA(λ) is the corresponding scaling function with the condition of SASA(1) = 1. This SASA term plays two inter-related roles: (1) it can effectively vary the surface tension surrounding the solute molecule and so indirectly perturbs near-solute solvent−solvent interactions; (2) correspondingly in the response order parameter Fλ, there is an additional (dSASA(λ)/dλ)Usasa(XN) term that allows nearsolute solvent−solvent fluctuations to be indirectly accelerated through the OSS treatment. More importantly, we generalize the response order parameter Fλ [Uss(XN) + (dSE(λ)/dλ)Use(XN) + (dSASA(λ)/

∂Hϕ1, ϕ2

(ϕ1,ϕ2)

∮[ϕ

m

∂Φ

dr λ′

(10)

in which ⟨∂Hϕ1,ϕ2/∂Φ⟩λ′ denotes (⟨∂Hϕ1,ϕ2/∂ϕ1⟩λ′, ⟨∂Hϕ1,ϕ2/ ∂ϕ2⟩λ′), a (ϕ1,ϕ2) space vector, and r represents a possible (ϕ1,ϕ2) space pathway between [ϕ1,min(λ′), ϕ2,min(λ′)] and (ϕ1,ϕ2). Following the strategy in a recent ABF implementation,35 the Monte Carlo integration method was used to ϕ1,ϕ2 calculate ∮ [ϕ ⟨∂Hϕ1,ϕ2/∂Φ⟩λ′ dr; this allows us to 1,min(λ′),ϕ2,min(λ′)] take into account all of the possible pathways in a maximum likelihood manner. Generalizing eq 5, we have the following equation to calculate the calibration constant Cgm for each λ′ state: 43


Article


Cg (λ′) = −max[gm(λ′, ϕ1 , ϕ2)] m

[ϕ1,max (λ′), ϕ2,max (λ′)]

∬[ϕ

− kTo ln

1,min(λ′), ϕ2,min(λ′)]

⎛ g (λ′, ϕ , ϕ ) − max[g (λ′, ϕ , ϕ )] ⎞ 1 2 1 2 m ⎟⎟ dϕ1 dϕ2 exp⎜⎜ m kT ⎝ ⎠ o

Then the target three-dimension function gm(λ,ϕ1,ϕ2) can be directly assembled from calculated two-dimension gm(λ′,ϕ1,ϕ2) functions. Generalizing eq 6, we have the following equation to adaptively calculate f m(λ):

used to treat long-range electrostatic interactions; in real space energy and force calculations, short-range electrostatic interactions were switched off at 12 Å. All the simulations were performed under the NPT ensemble. The temperature was set as 300 K, and the pressure was set as 1 atm. The Nóse-Hoover thermostat44 was utilized to maintain the constant temperature, and the Langevin piston algorithm45 was employed to maintain the constant pressure. The SHAKE algorithm46 was used to constrain all the bonds that involve water hydrogen atoms. The MD simulation time step was set as 1 fs.

⎛ ϕ1,max (λ′) ϕ2,max (λ′) gm(λ′, ϕ1, ϕ2) dϕ2 dϕ1 ⎜ ∫ϕ (λ′) ϕ1 ∫ϕ (λ′) exp kTo 2,min ⎜ 1,min − ϕ (λ′) ϕ (λ′) g (λ′, ϕ , ϕ ) λmin ⎜ ⎜ ∫ 1,max ∫ 2,max exp m kT 1 2 dϕ2 dϕ1 o ⎝ ϕ1,min(λ′) ϕ2,min(λ′)

∫

λ

ϕ

+

(λ′)

ϕ

( ( exp( exp(

(λ′)

∫ϕ 2,max(λ′) ϕ2 ∫ϕ 1,max(λ′) 2,min

ϕ2,max (λ′)

∫ϕ

2,min(λ′)

1,min

ϕ1,max (λ′)

∫ϕ

1,min(λ′)

) ) ) dϕ dϕ ⎞⎟⎟ dλ′ ) dϕ dϕ ⎟⎟⎠

gm(λ′, ϕ1, ϕ2) kTo

gm(λ′, ϕ1, ϕ2) kTo

1

1

(11)

2

2

(12)

The same as eq 6, eq 12 is an approximate way of estimating f m(λ); however it is sufficient for the purpose of predictive sampling. Generalized Orthogonal Space Tempering Simulation Result Analysis. For any collective variable set S of interest, the free energy surface can be constructed based on the following equation, Go(S) = −kTo ln ∑ e Hm(t ) − Ho(t )/ kTo δ(S[XN (t )] − S) t

+ constant

(13)

where t represents scheduled sample-collection time-steps in the near-equilibrium phase of the simulation. The above extended-ensemble-to-canonical-ensemble reweighting strategy has been employed in our earlier OSS based simulation25,36 and on-the-path random walk based simulation studies.37 It is a generalized form of the umbrella sampling formula,38 which was originally derived in the canonical-ensemble-to-canonicalensemble reweighting framework.

Figure 1. Deca-alanine model. Both the N-terminus and the Cterminus are neutrally capped. It is shown as an α-helix structure, in which, totally, six helical hydrogen bonds exist.

■

Details on the Generalized Orthogonal Space Tempering Simulation. The gOST algorithm was implemented in our customized version of the CHARMM program.47,48 In this study, we selected all the intrapeptide torsional terms (including the CMAP terms) and all the electrostatic energy terms into Uss(XN) and all the electrostatic interactions between the peptide and the solvent environment (including water molecules and ions) into Use(XN); in the Results and Discussion section, these two energy components are respectively denoted as UPP and UPW. For SE(λ), we picked √λ for no particular reason other than that it is a commonly employed functional form.28−31 For Usasa(XN), we used a surface-accessible-surface-area energy form49 in the CHARMM program. We collected experimentally determined temperaturedependent water surface tension changes γ(Tm) and set SASA(λ) as a B-spline function connecting the collected-points based on the relationship SASA(λ) = γ(300K/λ2)/γ(300K). To constrain λ within our preset range {0.801, 1.001}, following the original OST method,27 a θ-dynamics strategy was used. Specifically, we explicitly propagate a θ particle that travels periodically between −π and π; here, θ is the variable of

COMPUTATIONAL DETAILS To test the gOST method, the explicitly solvated deca-alanine (Ala10) peptide (Figure 1) was chosen as the model system. Although Ala10 is one of the most commonly studied systems, as far as we know, before this work, there has been no report on any successful “predictive” simulation of repetitive unfolding/ refolding of the explicitly solvated Ala10, in particular through a single continuous trajectory. As a matter of the fact, based on the gOST simulation result, our analysis suggests that commonly employed “predictive” sampling methods may not be effective on this seemingly simple system. Model and General Molecular Dynamics Simulation Setups. In the model, we have a deca-alanine peptide chain (Figure 1) solvated in a truncated octahedral water box; potassium and chloride ions were added to have the ionic strength adjusted around 0.15 M. The deca-alanine peptide and ions were modeled with the CHARMM2239/CMAP40 force field, and water molecules are treated with the modified TIP3P model.41,42 The particle mesh Ewald (PME) method43 was 44


Article


Figure 2. (a) Time-dependent fluctuations of the essential solute−solute interaction energy. (b) The time-dependent fluctuations of the helicity deviation collective variable.

Figure 3. (a) The free energy surface along the helicity deviation. (b) The two-dimensional free energy surface along the helicity deviation and the essential solute−solute interaction energy.

the λm(θ) function: λm(θ) = λmin + (λmax − λmin)λ(θ), in which λ(θ) takes the functional form as in the original OST method. Considering this θ-dynamics element, the expanded Hamiltonian in eq 7 should be

which was independently perturbed in the expanded Hamiltonian, fluctuates across several distinct regions. Each of these regions is featured with the corresponding enduring dynamics and undergoes possible transitions to its neighboring regions via rare fluctuations. To characterize the conformation corresponding to each of these UPP regions, based on our detailed survey of the trajectory, we defined a helicity deviation 2 1/2 collective variable, (Σi=6 i=1(di − 3.3) /6) , in which di represents the distance between a pair of α-helix hydrogen bond partners. As comparatively shown, the fluctuation of the helicity deviation collective variable (Figure 2b) is closely synchronous with the fluctuation of UPP (Figure 2a). Specifically, the lowest UPP region corresponds to the folded helical state with the lowest helicity deviations, the highest UPP region corresponds to the fully unfolded states with the highest helicity deviations, and the intermediate UPP regions correspond to different degrees of partially unfolded conformations. Further, the close correlation between the helicity deviation and the essential intrasolute interaction UPP can be directly revealed by the free energy surface along these two collective variables (Figure 3b). It should be noted that in this study, all the free energy surfaces were generated based on the reweighting formula of eq 13. In the gOST method, UPP, labeled as Fintra in eq 9, is employed as λ an independent response order parameter; the excellent correspondence between its fluctuations and essential conformational transitions demonstrates the advantage of the generalization treatment in the gOST formulation.

Hλ = λm(θ )Uss(XN) + SE(λm(θ ))Use(XN ) + [SASA(λm(θ )) − 1]Usasa(XN ) + Ue(XN ) N

+

∑ Pi 2 /2mi + pθ 2 /2mθ i=1

(14)

In our gOST simulation, the orthogonal space sampling temperature TES was set to 3000 K, corresponding to α = 0.9; the masses of the fictitious particles θ, ϕ1, and ϕ2 were respectively set as 2000, 0.0002, and 0.0002 atomic mass unit (amu); the restraint force constants Kϕ1 and Kϕ2 were set as 0.1 kcal −1 mol −1 . In gOST recursions, both f m (λ) and inter gm(λ,Fintra λ ,Fλ ) were updated every 1 ps (ps). The trajectories were saved every 1.0 ps for postprocessing analysis. For all the analysis based on eq 13, only the data after 60 ns (ns) were used; here the first 60 ns portion of the simulation serves to equilibrate the system and generate converged biasing energy functions.

■

RESULTS AND DISCUSSION The gOST Simulation Enabled Repetitive Ala 10 Unfolding and Refolding. As shown in Figure 2a, during the gOST simulation, UPP, the essential intrasolute energy, 45


Article


Figure 4. (a) Time dependent fluctuations of the end-to-end distance. (b) The two-dimensional free energy surface along the helicity deviation and the end-to-end distance.

method enabled successful exploration of a presumably more rugged energy landscape, we expect this sampling algorithm to be more effective when more accurate force fields are applied. The End-to-End Distance is A Poor Order Parameter. To describe helix−coil transitions, the end-to-end distance is the most commonly employed order parameter. As shown by earlier studies,51−53 fluctuations of this collective variable can closely follow Ala10 gas phase folding/unfolding transitions and thus be an effective sampling order parameter. In this predictive sampling based study, no system-specific collective variable was applied to guide the sampling. Therefore, the end-to-end distance can be evaluated in an unbiased manner. As comparatively shown in Figures 4a and 2b, in aqueous solution, the end-to-end distance does not monotonically change in correspondence to any folding or unfolding transition. Instead, independent of the conformational region being visited, it mostly fluctuates up and down around 14.8 Å, which is the endto-end distance of the fully folded helical structure. Interestingly, the fluctuation amplitude has some degree of correlation with the folding/unfolding transition. At the folded helical state, the fluctuation range is smaller, at the fully unfolded states, the fluctuation ranges are much larger, and at partially folded states, the fluctuation has medium magnitudes. This unique correlation behavior reflects the fact that unfolded regions of the peptide are very flexible and do not particularly prefer more extended conformations as they do in the gas phase. Apparently in aqueous solution, the peptide can access a much larger conformational space than in a vacuum. This can be further revealed by the two-dimensional free energy surface as the function of the helicity deviation and the end-to-end distance. Along the end-to-end distance direction, the free energy profile displays a largely pseudosymmetric shape with 14.8 Å as the symmetry center. At the states of “N-1” and “N2,” the disordered N-terminus fragments are by and large evenly distributed between extended and compact configurations. As an exception, when the partial unfolding occurs in the C-terminus end, the end-to-end distance has narrower and more specific distributions: at the C-1 state, the unfolded end is more stretched, and at the C-2 state, the unfolded end prefers to interact with the folded part of the peptide and so the region with a shorter end-to-end distance is more populated. When the helicity deviation value is larger than 4.0 Å, Ala10 is fully unfolded; interestingly, the peptide chain favors more compact coil conformations rather than stretched coil conformations. When fully unfolded, the conformation distribution is mainly

Further based on eq 13, we constructed the one-dimension free energy surface along the helicity deviation collective variable (Figure 3a). The free energy surface displays multiple minima with the global one corresponding to the native αhelical state. In the region with the helicity deviation lower than 4.0 Å, each of the local minima can be identified to correspond to a partially unfolded state. For instance, around 3.1 Å, the second lowest free energy minimum represents the state that loses the two N-terminus helical hydrogen bonds; here, this state is labeled as “N-2.” In regions with larger helicity deviation values, the conformational states are likely degenerate. Therefore, it is difficult to assign individual states directly from the one-dimension free energy surface. Nevertheless, overall, the helicity deviation is an excellent order parameter to describe Ala10 folding/unfolding transitions. We would like to note that in the following, we label various conformational states using symbols such as “N-l,” “C-m,” or “N-l/C-m;” taking “N-l/C-m” as the example, it denotes the partially unfolded state that loses l numbers of hydrogen bonds at the N-terminus end and m numbers of hydrogen bonds at the C-terminus end. Here, we define 3.45 Å as the cutoff value for the forming and breaking of an α-helix hydrogen bond. In our model, the CHARMM22/CMAP potential is used as the energy function. As has been known, this energy function likely overstabilizes the α-helical state and cannot accurately model generally expected cooperative helix−coil transitions.50 Our result agrees well with these general understanding. Specifically, the result shows that there is no clear broad free energy well around the unfolded state; instead, the CHARMM22/CMAP force field favors the fully folded state and also some of the partially unfolded states, such as the “N-2” state. Therefore, along folding/unfolding pathways, the peptide is very likely to be locally trapped by some of these over stabilized intermediate states; consequently, stepwise rather than concerted folding/unfolding pathways are more dominant. Due to its deficiency, this commonly used protein force field is not ideal for modeling peptide helix−coil transitions. However, because CHARMM22/CMAP features a more rugged energy surface than the true energy surface, the Ala10 model imposes a greater sampling challenge, and it can arguably serve as a better model system to test sampling techniques. As far as we know, before this work, there has been no report on any successful simulation of repetitive unfolding/refolding of the explicitly solvated Ala10, in particular through a single continuous trajectory. Based on the fact that in this study, the gOST 46


Article


Figure 5. (a) The two-dimensional free energy surface along the helicity deviation and the concertedness collective variable. (b) The thermodynamic network. “Native” labels the fully folded α-helix state.

Figure 6. (a) The two-dimensional free energy surface along the essential solute−solvent interaction energy and the essential solute−solute interaction energy. (b) The two-dimensional free energy surface along the helicity deviation and the solvent-accessible surface area dependent energy.

located between the region around ∼15.0 Å of the end-to-end distance and the region around ∼5.0 Å of the end-to-end distance. All the previous geometry-collective-variable-based MD studies had their sampling focused in the region with larger end-to-end distances, which correspond to extended coil conformations. Such discrepancy reveals the importance of the “predictive” sampling strategy, which does not rely on any preassumed guess on target systems. Obviously in aqueous solution, the end-to-end distance is a poor collective variable to describe Ala10 helix−coil transitions. The same conclusion was also drawn in a recent study based on the CHARMM36 potential energy function.53 Again, our analysis was based on a “predictive” sampling simulation, which, without the guidance of any system-specific order parameter, explores a much broader range of the end-to-end distance and thus provides more accurate and unbiased information for the assessment analysis. The CHARMM22/CMAP Energy Function Leads to Stable Partially-Unfolded Intermediates. As mentioned earlier, the CHARMM22/CMAP energy potential biases the explicitly solvated Ala10 to favor clearly stepwise rather than commonly expected concerted folding/unfolding pathways. Seeking to understand the cause, we define an order parameter, “concertedness,” to differentiate concerted and stepwise pathways. Specifically, concertedness is defined as the root-

2 1/2 mean-square deviation, (∑i=6 i=1(di − d̅) /6) , in which d̅ stands for the average value of the distances between each pair of the α-helix hydrogen bond partners. Here, a lower concertedness value corresponds to a structure along a more concerted pathway and a higher value corresponds to a structure along a more stepwise pathway. As shown by the two-dimensional free energy profile with the helicity deviation and the concertedness as the collective variables (Figure 5a), for any plausible transition between the fully folded helical state and the completely unfolded state, the pathway has to go through certain partially unfolded states with relatively high “concertedness” values, i.e. conformational intermediates that characteristically represent stepwise transitions, such as states “N-3” and “C-3” etc. This is due to the fact that these intermediates and their precursor intermediates, such as the “N-2” and “C-2” states, are too stable. To more comprehensively view folding/unfolding pathways, based on the obtained samples, we constructed the thermodynamic network view of the free energy landscape (Figure 5b). In this network landscape, the area of each sphere is proportional to the population of the corresponding state and the width of each line is proportional to the number of the transitions between each pair of the states. The thermodynamic network reveals that there is no direct connection between the fully folded helical state and the completely unfolded (N-3/C-

47


Article


Figure 7. (a) Time dependent changes of surface-proximity-dependent water number profiles. (b) Time dependent changes of the number of waters in the layer of 1.25 Å above the molecular surface. (c) Time dependent changes of the number of waters in the layer of 2.05 Å above the molecular surface. (d) Time dependent changes of the number of waters in the layer of 5.05 Å above the molecular surface. The green lines in a, b, and c represent the local averages of the number of water.

process, the essential solute−solvent and solute−solute interactions always intimately interplay. Specifically, the solute−solvent and solute−solute interaction energies not only fluctuate in the opposite direction but also closely compensate each other. This observation suggests that solute conformational changes require synchronous solvent responses, which are featured by the compensating fluctuations of essential solute−solvent and solute−solution interactions. Such observation further confirms the major motivation underlying the gOST algorithm design: separately treating the solute−solute and solute−solvent interaction components in the response order parameter set. In addition to solute−solvent interactions, solvent−solvent interactions are also likely involved in solute conformational transitions. In gOST, a solvent-accessiblesurface area (SASA) dependent term is introduced to accelerate near-solute solvent response fluctuations. As shown by Figure 6b, the SASA energy and the helicity deviation are also highly correlated; such close correspondence implies the importance of the SASA term in the gOST formulation. Solvent Cooperative Transitions Play a Pivotal Role in Ala10 Folding/Unfolding Processes. To understand how solvent structural changes couple with Ala10 conformational transitions, we calculated the surface-proximity-dependent water number profile for each collected sample and monitored its time-dependent changes. The molecular surface definition and the surface-proximity-dependent water number calculation are based on the work of Willard and Chandler.54 In this work, each water layer is defined as a shell with the thickness of 0.1 Å. As shown in Figure 7a, in correspondence with solute structure

3) state. Interestingly, the completely unfolded (N-3/C-3) state is mainly connected with the “N-3” and “C-3” states; this suggests that only for a short fragment with less than three helical hydrogen bonds could the concerted folding/unfolding transition be a dominant pathway. In order for the solvated Ala10 peptide at either the “N-3” or “C-3” state to become fully folded, it has to follow a stepwise path, respectively going through the “N-2” state or “C-2” state, both of which can form local thermodynamic traps. Moreover, our result shows that the C-terminus end is more stable than the N-terminus end. When the unfolding process is initiated from the N-terminus end, the system is most likely to be trapped at the overstabilized N-2 and its accessory states, such as the N-2/C-1 and N-2/C-2 states, before proceeding to the fully unfolded state. The above analysis is expected to provide insightful details and guidance for future force field reparameterization and refinement efforts. For the fact that the gOST algorithm can enable detailed sampling of Ala10 folding/unfolding processes, this method can serve as an essential part of a future routine procedure for protein energy function assessment. Compensating Fluctuations of Solute−Solvent and Solute−Solute Interactions Drive Folding/Unfolding Transitions. As discussed earlier, fluctuations of the essential solute−solute interaction energy Upp are highly correlated with Ala10 conformational changes (Figure 3b). In order to understand how solvent molecules influence helix−coil transitions, we constructed the two-dimensional free energy surface along both the essential solute−solvent interaction energy UPW and Upp (Figure 6a). Notably throughout the whole 48


Article


Figure 8. (a) Time dependent fluctuations of the total potential energy. (b) The two-dimensional free energy surface along the helicity deviation and the total potential energy. (c) The two-dimensional free energy surface along the helicity deviation and the total solute torsion energy. (d) The twodimensional free energy surface along the helicity deviation and Upp + (1/2)√λUPW.

method strongly depends on how well fluctuations of a prechosen energy order parameter correlate with the progress of molecular events. With a poor correlation, projected fluctuation enlargements on essential degrees of freedom are expected to be small, and sampling efficiency is expected to be low. As discussed above, fluctuations of both of the response order parameters in gOST strongly correlate with unfolding/ folding transitions (Figures 3b and 6a). This should naturally explain the effectiveness of the gOST method on the sampling of the solvated Ala10 model system. Figure 8a shows that throughout the simulation, during which folding and unfolding events repetitively occur, the total potential energy does not display any significant change. As discussed earlier, during solute conformational transitions, distinct interaction components do not always fluctuate in the same direction. Instead, in this study, the essential solute− solute and solute−solvent interactions display a compensating fluctuation behavior; consequently, the fluctuation magnitude of the overall potential energy is very moderate. The twodimensional free energy surface along the helicity deviation and the total potential energy (Figure 8b) reveals that there is no obvious correlation between conformational changes and total potential energy fluctuations. Due to its poor correlation with the conformational change, when the total potential energy is activated, only a small fraction of the fluctuation enlargement can be effectively projected onto the degrees of freedom essential to the conformational changes. Thus, methods relying on the total potential energy fluctuation promotion, such as TREM, ST, and the original aMD method, are unlikely to be effective on Ala10 folding/unfolding transitions. As an alternative, in some predictive sampling practices, the total solute torsional energy (corresponding to the sum of the

changes, water structures undergo characteristic collective transitions, for instance, in the layers that are 1.25 and 5.05 Å above the molecular surface. Taking a close look at the solvent layer that is 1.25 Å above the molecular surface (Figure 7b), we can tell that instead of slowly relaxing to equilibrium configurations after major solute conformational changes, water molecules undergo simultaneous configuration changes Such simultaneous transitions also occur in more distant layers, for instance the one that is 5.05 Å above the molecular surface (Figure 7d). Though with smaller amplitudes, the structural fluctuations of the distant layers, particularly from the local averaging view, display almost identical features as the fluctuations in the 1.25 Å layer. Even for the layers that undergo very small structural changes, for instance the one that is 2.05 Å above the molecular surface (Figure 7c), a similar transition pattern is displayed as well. The above analysis reveals that near-solute solvent molecules are an integral part of the solute system, and water cooperative fluctuations occur simultaneously with and are essential to slow solute conformational transitions. Therefore, in a sampling algorithm design, it is dangerous to assume that solvent structures can always quickly relax and so empirically designate the solute region as the only “hot” area for active sampling treatment and the solvent region as the “cold” area to passively respond. Fluctuation Behaviors of Common Potential Energy Based Order Parameters. As discussed in the Introduction, common “predictive” sampling methods, including TREM, HREM, ST, SS, and accelerated MD (aMD)55 etc. were formulated to enlarge fluctuations of certain potential-energy based order parameters with the hope that the promoted interaction fluctuations can facilitate rare event barrier crossings. In practice, the effectiveness of a “predictive” sampling 49


Article

Journal of Chemical Theory and Computation total torsional energy and the total CMAP energy in this study) is employed as the sampling order parameter. Figure 8c shows that there is no obvious correlation between total torsional energy fluctuations and helix−coil transitions. Therefore, we expect that enlarging the fluctuation of the total solute torsional energy alone is unlikely to enable effective sampling of the solvated Ala10 peptide. It is worth noting that in our expanded Hamiltonian treatment, the total torsional energy is included in Uss(XN). From the physical viewpoint, conformational transitions are governed by the fluctuations of both nonbond interaction energies and torsional energies. Together with the discussion in the previous subsection, we can tell that Ala10 folding/unfolding transitions are mainly driven by solute− solute and solute−solvent electrostatic interaction fluctuations. The HREM and ST methods are based on the general expanded Hamiltonian as described by eq 1, which is employed as the basis of the perturbation treatment in this study. Specifically, we chose λ1/2 for SE(λ). Here, we evaluate Upp + (1/2)λ1/2UPW, the fluctuation of which is to be enlarged by the HREM or ST method basesd on eq 1. As shown by the free energy surface in Figure 8d, fluctuations of Upp + (1/2)λ1/2UPW do not correlate well with changes of the helicity deviation. This particular expanded Hamiltonian (eq 1) has been increasingly popular largely due to the fact that it can produce relatively high acceptance ratios and therefore high λ-space fluxes.31 Our study shows that essential conformational changes may not be sensitive to the fluctuation of Upp + (1/2)λ1/2UPW; then due to a poor coupling with the configuration space response, even a high λ-space flux cannot lead to effective sampling of solute conformational transitions. Therefore, such methods are unlikely to be effective on Ala10 folding/unfolding transitions.

importantly, our analysis suggests that solvent cooperative fluctuations intimately interplay with Ala10 conformational transitions. For the fact that folding/unfolding transitions are mainly driven by the compensating fluctuations of essential solute−solvent and solute−solute interactions, as indicated by our assessment analysis result, commonly employed “predictive” sampling methods may not be effective on this seemingly “simple” system. Although this paper is to present a new sampling method, we would also like to share our insight regarding how to sample microscopic events that are driven by compensating fluctuations of distinct energy components. If a certain sum of these components, such as total potential energy or total interaction (both intramolecular and intermolecular interaction) energy related to the solute molecule, is employed as a sampling order parameter, expected large degeneracy can lead to the so-called “entropy-barrier” sampling issue, which may introduce large diffusion sampling overhead and greatly lower the sampling effectiveness. In the orthogonal space sampling scheme, we can separately treat their compensating components in the response order parameter set and enable robust activation of configuration responses. Finally, we hope to use this design and demonstration of the gOST method to illustrate the power the OSS scheme in physics-based sampling method developments.

CONCLUDING REMARKS In aqueous solution, solute conformational transitions are governed by intimate interplays of the fluctuations of solute− solute, solute−water, and water−water interactions. To predictively sample solute conformational changes in aqueous solution, a common strategy is to construct an expanded Hamiltonian through a series of Hamiltonian perturbations and thereby broaden the distributions of certain interactions of focus. As is increasingly realized, general expanded Hamiltonian methods suffer from the issue of slow configuration response, particularly when solute structural changes are intimately coupled with solvent rearrangements. To more effectively tackle this issue, in this work, following the orthogonal space sampling scheme, we designed a generalized orthogonal space tempering (gOST) method. In this sampling method design, the emphasis is to separately treat the solute−solute and solute−solvent components in the response order parameter set and implicitly perturb near-solute water−water fluctuations. To demonstrate the gOST method, we carried out a molecular dynamics simulation study on the explicitly solvated deca-alanine (Ala10) peptide. Based on a fully automated sampling protocol, the gOST simulation enabled repetitive folding and unfolding of the solvated peptide within a single continuous trajectory and allowed for detailed constructions of Ala10 folding/unfolding free energy surfaces. The gOST result allows us to identify the helicity deviation collective variable that can reasonably describe the Ala10 helix−coil transition. In addition, our result shows that the CHARMM22/CMAP energy function overstabilizes certain partially unfolded states that lead to stepwise folding/unfolding processes. More

Present Address

■

AUTHOR INFORMATION

Corresponding Authors

*Tel.: 850-645-6884. Fax: 850-644-7244. E-mail: lzheng@fsu. edu. *Tel.: 850-645-6884. Fax: 850-644-7244. E-mail: yyang2@fsu. edu.

■

∥

Washington University in St. Louis, Campus Box 1134, 1 Brookings Drive, St. Louis, MO 63130-4899 Author Contributions §

Equal contribution to this work

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Funding supports from the National Science Foundation (MCB1158284) and the National Institutes of Health (R01GM111886) are acknowledged. We would also like to thank the Oak Ridge National Laboratory Supercomputing Center and the Florida State University Research Computing Center for the computing support.

■

REFERENCES

(1) Berne, B. J.; Straub, J. E. Novel methods of sampling phase space in the simulation of biological systems Curr. Curr. Opin. Struct. Biol. 1997, 7, 181−189. (2) Okamoto, Y. Generalized-ensemble algorithms: Enhanced sampling techniques for Monte Carlo and molecular dynamics simulations. J. Mol. Graphics Modell. 2004, 22, 425−439. (3) Elber, R. Long-timescale simulation methods. Curr. Opin. Struct. Biol. 2005, 15, 151−156. (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) Yang, W.; Nymeryer, H.; Zhou, H. X.; Berg, B. A.; Brüschweiler, R. Quantitative computer simulations of biomolecules: A snapshot. J. Comput. Chem. 2008, 29, 668−672. 50


Article

Journal of Chemical Theory and Computation (6) Zuckerman, D. M. Equilibrium sampling in biomolecular simulation. Annu. Rev. Biophys. 2011, 40, 41−62. (7) Abrams, C.; Bussi, G. Enhanced sampling in molecular dynamic using metadynamics, replica-exchange, and temperature-acceleration. Entropy 2014, 16, 163−199. (8) Swendsen, R. H.; Wang, J. S. Replica Monte Carlo simulation of spin-glasses. Phys. Rev. Lett. 1986, 57, 2607−2609. (9) Geyer, C. J. Practical Markov chain Monte Carlo. Stat. Sci. 1992, 7, 473−483. (10) Hukushima, K.; Nemoto, K. Exchange Monte Carlo method and application to spin glass simulations. J. Phys. Soc. Jpn. 1996, 65, 1604− 1608. (11) Hansmann, U. H. E. Parallel tempering algorithm for conformational studies of biological molecules. Chem. Phys. Lett. 1997, 281, 140−150. (12) Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 1999, 314, 141−151. (13) Sugita, Y.; Kitao, A.; Okamoto, Y. Multidimensional replicaexchange method for free-energy calculations. J. Chem. Phys. 2000, 113, 6042−6051. (14) 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. (15) 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. (16) Liu, P.; Huang, X.; Zhou, R.; Berne, B. J. Hydrophobic aided replica exchange: an efficient algorithm for protein folding in explicit solvent. J. Phys. Chem. B 2006, 110, 19018−19022. (17) 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−1783. (18) Marinari, E.; Parisi, G. Simulated tempering: A new Monte Carlo scheme. Europhys. Lett. 1992, 19, 451−458. (19) Irback, A.; Potthast, F. Studies of an off-lattice model for protein folding: Sequence dependence and improved sampling at finite temperature. J. Chem. Phys. 1995, 103, 10298−10305. (20) Li, H.; 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. (21) Duane, S.; Kennedy, A. D.; Pendleton, B. J.; Roweth, D. Hybrid Monte Carlo. Phys. Lett. B 1987, 195, 216−222. (22) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182−7190. (23) Nosé, S.; Klein, M. L. Constant pressure molecular dynamic for molecular systems. Mol. Phys. 1983, 50, 1055−1076. (24) Kong, X.; Brooks, C. L., III λ-dynamics: A new approach to free energy calculations. J. Chem. Phys. 1996, 105, 2414−2423. (25) Zheng, L. Q.; Chen, M. E.; Yang, W. Simultaneous escaping of explicit and hidden free energy barriers: Application of the orthogonal space random walk strategy in generalized ensemble based conformational sampling. J. Chem. Phys. 2009, 130, 234105. (26) Zheng, L. Q.; Chen, M. E.; 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. (27) Zheng, L. Q.; Yang, W. Practically Efficient and Robust Free Energy Calculations: Double-Integration Orthogonal Space Tempering. J. Chem. Theory Comput. 2012, 8, 810−823. (28) Su, L.; Cukier, R. I. Hamiltonian replica exchange method study of Escherichia coli and Yersinia pestis HPPK. J. Phys. Chem. B 2009, 113, 16197−16208. (29) Moors, S. L. C.; Michielssens, S.; Ceulemans, A. Improved replica exchange method for native-state protein sampling. J. Chem. Theory Comput. 2011, 7, 231−237.

(30) Terakawa, T.; Kameda, T.; Takada, S. On easy implementation of a variant of the replica exchange with solute tempering in GROMACS. J. Comput. Chem. 2011, 32, 1228−1234. (31) Wang, L.; Friesner, R. A.; Berne, B. J. Replica exchange with solute scaling: A more efficient version of replica exchange with solute tempering (REST2). J. Phys. Chem. B 2011, 115, 9431−9438. (32) Lv, C.; Zheng, L. Q.; Yang, W. Generalized essential energy space random walks to more effectively accelerate solute sampling in aqueous environment. J. Chem. Phys. 2012, 136, 044103. (33) Darve, E.; Pohorille, A. Calculating free energies using average force. J. Chem. Phys. 2001, 115, 9169−9183. (34) Sitkoff, D.; Sharp, K.; Honig, B. Accurate calculation of hydration free energies using macroscopic solvent models. J. Phys. Chem. 1994, 98, 1978−1988. (35) Hénin, 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. (36) Min, D.; Chen, M. E.; Zheng, L. Q.; Jin, Y.; Schwartz, M. A.; Sang, Q. X. A.; Yang, W. Enhancing QM/MM Molecular Dynamics sampling in explicit environments via an orthogonal-space-randomwalk-based strategy. J. Phys. Chem. B 2011, 115, 3924−3935. (37) Cao, L.; Lv, C.; Yang, W. Hidden conformation events in DNA base extrusions: A generalized ensemble path optimization and equilibrium simulation study. J. Chem. Theory Comput. 2013, 9, 3756−3768. (38) Torrie, G. M.; Valleau, J. P. Non-physical sampling distributions in Monte-Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23, 187−199. (39) MacKerell, A. D., Jr; Bashford, D.; Bellot, M.; Dunbrack, R. L., Jr; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCathy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (40) MacKerell, A. D., Jr; Feig, M.; Brooks, C. L., III Extending the treatment of backbone energetics in protein force fields: Limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J. Comput. Chem. 2004, 25, 1400−1415. (41) 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. (42) Neria, E.; Fischer, S.; Karplus, M. Simulation of activation free energies in molecular systems. J. Chem. Phys. 1996, 105, 1902−1921. (43) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An nLog(n) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089−10092. (44) Nóse, S. A unified formulation of the constant temperature molecular-dynamics methods. J. Chem. Phys. 1984, 81, 511−519. (45) Feller, S. E.; Zhang, Y. H.; Pastor, R. W.; Brooks, B. R. Constant-pressure molecular-dynamics simulation: The Langevin piston method. J. Chem. Phys. 1995, 103, 4613−4621. (46) Ryckaert, J.; 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. (47) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 1983, 4, 187−217. (48) Brooks, B. R.; Brooks, C. L., III; Mackerell, A. D., Jr; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Calfisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; 51


Article

Journal of Chemical Theory and Computation York, D. M.; Karplus, M. CHARMM: The biomolecular simulation program. J. Comput. Chem. 2009, 30, 1545−1614. (49) Lee, M. S.; Feig, M.; Salsbury, F. R., Jr; Brooks, C. L., III New analytic approximation to the standard molecular volume definition and its application to generalized Born calculations. J. Comput. Chem. 2003, 24, 1348−1356. (50) Best, R. B.; Mittal, J.; Feig, M.; MacKerell, A. D. Inclusion of many-body effects in the additive CHARMM protein CMAP potential results in enhanced cooperativity of α-Helix and β-Hairpin Formation. Biophys. J. 2012, 103, 1045−1051. (51) 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, 3559−3566. (52) Hénin, J.; Chipot, C. Overcoming free energy barriers using unconstrained molecular dynamics simulations. J. Chem. Phys. 2004, 121, 2904−2914. (53) Hazel, A.; Chipot, C.; Gumbart, J. C. Thermodynamics of decaalanine folding in water. J. Chem. Theory Comput. 2014, 10, 2836− 2844. (54) Willard, A. P.; Chandler, D. Instantaneous liquid interfaces. J. Phys. Chem. B 2010, 114, 1954−1958. (55) Hamelberg, D.; Mongan, J.; McCammon, J. A. Accelerated molecular dynamics: A promising and efficient simulation method for biomolecules. J. Chem. Phys. 2004, 120, 11919−11929.

52


Predictive Sampling of Rare Conformational Events in Aqueous

Recommend Documents