Zooming Across the Free-Energy Landscape, Shaving Barriers and

Subscriber access provided by CAL STATE UNIV SAN FRANCISCO

Biophysical Chemistry, Biomolecules, and Biomaterials; Surfactants and Membranes

Zooming Across the Free-Energy Landscape, Shaving Barriers and Flooding Valleys Haohao Fu, Hong Zhang, Haochuan Chen, Xueguang Shao, Christophe Chipot, and Wensheng Cai J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.8b01994 • Publication Date (Web): 03 Aug 2018 Downloaded from http://pubs.acs.org on August 4, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

Zooming across the Free-Energy Landscape, Shaving Barriers and Flooding Valleys Haohao Fu,† Hong Zhang,† Haochuan Chen,† Xueguang Shao,†,#,‡ Christophe Chipot,*,§,∥,⊥ and Wensheng Cai,*,†,# †

Research Center for Analytical Sciences, College of Chemistry, Tianjin Key

Laboratory of Biosensing and Molecular Recognition, Nankai University, Tianjin 300071, China #

Collaborative Innovation Center of Chemical Science and Engineering (Tianjin),

Tianjin 300071, China ‡

State Key Laboratory of Medicinal Chemical Biology, Tianjin 300071, China

§

Laboratoire

International

Associé

CNRS

and

University

of

Illinois

at

Urbana−Champaign, Vandœuvre-lès-Nancy F-54506, France ∥LPCT,

UMR 7019 Université de Lorraine CNRS, Vandœuvre-lès-Nancy F-54500,

France ⊥Department

of Physics, University of Illinois at Urbana−Champaign, 1110 West

Green Street, Urbana, Illinois 61801, United States Corresponding Author *Email: [email protected] (W. C.), [email protected] (C. C.) 1

ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 27

ABSTRACT: A robust importance-sampling algorithm for mapping free-energy surfaces over geometrical variables, coined meta-eABF, is introduced. This algorithm shaves

the

free-energy

barriers

and

floods

valleys

by

incorporating

a

history-dependent potential term in the extended adaptive biasing force (eABF) framework. Numerical applications on both toy models and nontrivial examples indicate that meta-eABF explores the free-energy surface significantly faster than either eABF or metadynamics (MtD) alone, without the need of stratifying the reaction pathway. In some favorable cases, meta-eABF can be as much as five times faster than other importance-sampling algorithms. Many of the shortcomings inherent to eABF and MtD, like kinetic trapping in regions of configurational space already adequately sampled, the requirement of prior knowledge of the free-energy landscape to set up the simulation, are readily eliminated in meta-eABF. Meta-eABF, therefore, represents an appealing solution for a broad range of applications, especially when both eABF and MtD fail to achieve the desired result.

TOC GRAPHICS

2


Page 3 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Enhanced sampling has been widely used to explore physicochemical processes, the timescale of which goes beyond the scope of brute-force molecular dynamics (MD).1-3 Generalized-ensemble4 and importance-sampling schemes5 form two classes enhanced-sampling algorithms. Generalized-ensemble simulations achieve random walk in configurational space by introducing an auxiliary probability weight factor. Conversely, in importance-sampling algorithms, sampling along some geometric variable is artificially encouraged, because it is believed to describe the process of interest. Umbrella sampling (US),6,7 metadynamics (MtD)8-11 and adaptive biasing force (ABF)12-14 are commonly employed importance-sampling algorithms. Each of them has advantages and drawbacks or limitations. For example, the original US algorithm6 requires prior knowledge of the free-energy profile, or potential of mean force (PMF), underlying the physicochemical process to define the umbrella potential, while in its common acceptation today,7 it stratifies the transition coordinate space into very small ranges, or windows, wherein the free-energy barrier can be overcome by thermal motion. Turning to overly small windows may, however, hamper sampling in the orthogonal space15,16 — a portion of configurational space that cannot be properly described by the chosen transition coordinate. MtD, an heir of the local elevation10 and conformational flooding11 algorithms, adopts a time-dependent potential composed of Gaussians to encourage exploration of the transition coordinate in regions of configurational space that have not been

3



sampled. Using this strategy, one can use larger windows, compared with US, and still expect suitable sampling efficiency. As can be easily imagined, however, the width and height of the Gaussians greatly affects the performance of the MtD algorithm in practice. Moreover, since MtD records no information about the orthogonal space, convergence of the simulations may be poor when the transition coordinate is strongly coupled to slow degrees of freedom in the orthogonal space.

In stark contrast, ABF flattens the free-energy landscape by adding biasing forces along the transition coordinate. In an ABF simulation, windows of a size comparable to that in MtD can be utilized, and no prior knowledge of the PMF is needed. However, the choice of the coarse variables must meet some prerequisites for the classical ABF algorithm.17 In addition, since estimating the biasing forces in poorly sampled areas of configurational space can prove difficult, exploration of the transition coordinate may get stuck in oversampled areas as a result of incorrect biases being added.

There are many variants of the aforementioned approaches. Usually, shortcomings of an algorithm can be partially overcome by its variants. For example, extended ABF (eABF) eliminates the prerequisites of classical ABF by introducing extended degrees of freedom, and improves the convergence rate of the simulation through coarse-grained dynamics.17,18 Extended-generalized ABF uses N one-dimensional free-energy profiles determined by eABF to approximate the N-dimensional free-energy landscape, thereby significantly reducing the computational cost of PMF

4


Page 4 of 27

Page 5 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


calculations.19 Multiple-walker ABF combines ABF with a replica-exchange scheme, in an attempt to achieve ergodic sampling in orthogonal space.20 Conversely, well-tempered MtD adds Gaussians that decrease in height over time, thus, making MtD converge smoothly.21,22 Bias-exchange MtD introduces a replica-exchange scheme into the MtD framework to improve the sampling efficiency, especially in orthogonal

space.23,24

Unified

free-energy

dynamics

combines

MtD

with

temperature-accelerated molecular dynamics to improve the convergence rate in high-dimensional free-energy calculations.25,26 Adaptively biased molecular dynamics reduces the numerical cost of MtD from O(t2) to O(t) by adding non-Gaussian terms.27 However, the broad range of importance-sampling algorithms available can be confusing to the layman because comparing objectively their merits and drawbacks, as well as their respective performance has often proven difficult in practice.

Different schools of thought are accustomed to different approaches; they seldom depart from their favorite algorithm, especially for nontrivial applications. No algorithm, however, has proven yet able to perform universally well. Under certain circumstances, poor choice of the method can lead to unreliable results (see movies in the Supporting Information). It is, therefore, of utmost importance to design a “general” importance-sampling algorithm, germane for the broadest possible range of applications of physical, chemical and biological relevance. Toward this end, we put forth meta-eABF, in which the history-dependent memory kernel of MtD is introduced in the eABF framework, while an umbrella integration-based15 or CZAR estimator18 is used to evaluate rigorously the gradient of the free energy. In this 5



Page 6 of 27

scheme, the advantages of MtD and eABF are combined synergistically to overcome their respective shortcomings. The sampling efficiency and convergence rate are greatly improved compared with plain eABF and MtD, as detailed below. Meta-eABF may, therefore, be considered as a first option of choice for the exploration of the free-energy landscape that underlies complex physicochemical processes.

The theoretical underpinnings of meta-eABF can be summarized as follows. In conventional eABF, a fictitious particle is coupled to the coarse variable by means of a stiff spring. The applied bias on the particle relies on the average spring force:

Fbias (ξ ' ) = K (ξ ' − ξ )

Ξ'

= K (ξ ' − ξ

Ξ'

),

(1)

where K is the force constant of the spring, ξ', the extended coarse variable, ξ, the “real” coarse variable. Ξ' represents the fictitious particle-restraining ensemble. The free-energy change along ξ', namely the extended PMF, ΔA’, is an approximation of the real PMF, ΔA. In eABF, accurate evaluation of Fbias for a region of configurational space requires appropriate sampling, which, in turn, requires correct evaluation of Fbias. Exploring an undersampled region of configurational space may, therefore, become difficult sometimes, as shown in Fig. 1A. In meta-eABF, an MtD-like memory kernel is incorporated into the extended system alongside the eABF biasing force, thus, leading to:

Fbias (ξ ) = Fbias, eABF (ξ ) + Fbias, MtD (ξ ) = K (ξ − ξ '

'

'

'

Ξ'

)+

dU MtD (ξ ' , t ) dξ '

,

(2)

where UMtD(ξ',t) is the time-dependent MtD-like memory kernel 8. The extended PMF 6


Page 7 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


then writes, ' ' . ΔA' = ΔAeABF + ΔAMtD

(3)

In other words, the eABF algorithm samples the extended free-energy surface minus the sum of MtD Gaussian distributions, while the MtD algorithm explores the extended

PMF

minus

the

eABF

bias.

Unlike

extended-Lagrangian-based

metadynamics (eMtD)28, both Gaussian- and mean-force-based biases, contributed respectively by MtD and eABF, are added in a meta-eABF simulation. Here, MtD-like Gaussians are incorporated using the same criterion as that in a conventional eMtD simulation, while the eABF mean force is evaluated on the free-energy surface as equal to ΔA’ - ΔA’MtD, where ΔA’ is the original extended free-energy and ΔA’MtD is derived from the MtD memory kernel. The method presented in our manuscript may, therefore, be described as “MtD plus eABF with an estimator” (meta-eABF). Although both ΔA’eABF and ΔA’MtD change as a function of time, the sum of the two remains a constant if the simulation has converged (see Figure S1 in the Supporting Information). Inheriting their respective merits, this combination of eABF and MtD can shave free-energy barriers and flood valleys at the same time, thereby obviating the need of an accurate evaluation of Fbias, since the memory kernel of MtD encourages sampling of the undersampled space, as shown in Figure 1C. The true PMF, ΔA, can be recovered in a meta-eABF simulation by means of an umbrella integration-based15 or CZAR estimator18. Figure 1 illustrates the basic principle of eABF, MtD and meta-eABF. The use of meta-eABF in NAMD29 is detailed in Figure

7



S2. With the recent implementation of an extended-Lagrangian-based ABF variant30,31 and the corresponding estimators17,18, all MD engines supporting the Colvars30 or Plumed32 free-energy modules can potentially run meta-eABF simulations. It ought to be noted that the fundamental principle of meta-eABF is different from that of the recently developed MtD-ABF hybrid method33 (see the Supporting Information for details).

Figure 1. Schematic illustration of eABF, MtD and meta-eABF. Through biasing the extended collective variables that are coupled to the real ones by means of harmonic 8


Page 8 of 27

Page 9 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


springs, eABF and meta-eABF explore a mollified free-energy surface. See ref. 17 and 18 for more details. As shown in Figure 1, combination of eABF and MtD results in a dual acceleration, which translates to the following advantages:

(a) In an eABF simulation, exploration of the undersampled regions of configurational space may be slow, on account of the inherent difficulty to estimate the biasing force, namely Fbias. In meta-eABF, however, the Gaussian potentials will push the system to these undersampled regions (see Figure 1C). Moreover, the MtD-like memory kernel renders sampling across the reaction pathway uniform, while the sampling distribution of other ABF-based algorithms squarely relies on the quality of the estimation of Fbias.

(b) The quality of an MtD simulation highly depends on the height and width of the Gaussians. If the height and width are small, the time needed to fill up the valleys can be exaggeratedly long. Conversely, if they are too large, the molecular assembly will be far from equilibrium at the initial stage of the simulation, and the latter may become unstable due to the magnitude of the biases, hence, lowering the reliability of the PMF accordingly. In meta-eABF simulations, the extended free-energy surface is partially flattened by the biasing force of eABF, as shown in Figure 1C. Since the free-energy barriers are lowered after the corresponding region is adequately sampled, the valleys that ought to be filled are shallow. A small Gaussian height and width, therefore, appears to constitute a safe choice. 9



(c) By using a small Gaussian height and width, meta-eABF can sample the free-energy landscape efficiently and recover rigorously the PMF by means of an umbrella integration-based15 or CZAR estimator18. Although MtD can use a large Gaussian height and width to achieve high sampling efficiency, the PMF cannot always be recovered with the desired accuracy. This issue is, however, addressed to some extent in well-tempered MtD.21 An additional parameter that relies on prior knowledge of the PMF, namely the bias temperature, is, nevertheless, needed with this algorithm.

(d) Large windows can be employed in meta-eABF, which has proven to enhance significantly orthogonal-space sampling. For instance, when two parallel valleys of the free-energy surface are separated by a high barrier in orthogonal space, if a strategy consisting of stratifying the reaction pathway in small windows is used, the least-free-energy pathway connecting the two local minima is likely to be interrupted by the boundaries of the windows, as shown in Figure S3.

To illustrate the sampling efficiency of meta-eABF, simulations characterizing the isomerization of N-acetyl-N’-methyl-L–alanylamide (NANMA) using eABF, MtD and meta-eABF were carried out. Traditionally, geometrical variables, specifically the backbone torsional angles φ and ψ have been used as the reaction-coordinate model,34-36 as shown in Figure 2A. Here, to investigate both the transition-coordinate space and orthogonal-space sampling efficiency, only torsional angle φ was chosen as the coarse variable, while the other backbone torsion, ψ, corresponds to the

10


Page 10 of 27

Page 11 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


orthogonal space. Sampling along (φ, ψ) was monitored, as shown in Figure 2C-F. The sampling by eABF covers the whole transition-coordinate space, but some of the local minima in orthogonal space are missing. MtD with a relatively small Gaussian height and width (Figure 2D) reveals an improved ability to sample the orthogonal space,

relative

to

plain

eABF.

However,

in

this

particular

instance,

transition-coordinate-space sampling by MtD is not as good as that of eABF. Either increasing the height and width of the Gaussians in MtD, or turning to meta-eABF can lead to a significant improvement in sampling efficiency, as shown in Figure 2E and 2F. Using the former scheme, however, may drive the simulation far from equilibrium and hamper convergence, as depicted by the realistic examples “shuttling of a rotaxane” and “translocation of a halide ion” detailed below. This shortcoming does not exist in meta-eABF simulations.

To clarify the idea of universality of meta-eABF, a simulation describing a more general case, namely a particle moving on a potential-energy surface constructed mathematically, was performed. The result clearly shows that, with a relatively small Gaussian height and width, meta-eABF has a sampling efficiency clearly superior to that of eABF and MtD (see Figure S4). The remarkable sampling efficiency of meta-eABF is, therefore, not specific to the NANMA toy model, but holds for the exploration of other, more complex free-energy surfaces. The parameters used in Figure 2F and Figure S4, namely hillWeight = 0.1 kcal/mol and hillWidth = 5 bin width in the Colvars vocabulary,30 appear to constitute a safe choice in meta-eABF simulations (See Figure S5). This set of parameters has proven appropriate for all the 11



examples reported in this work, and in a variety of unpublished applications examined by our respective research groups.

Figure 2. Simulation describing the isomerization of NANMA. (A) The backbone angles (φ, ψ). (B) The reference free-energy landscape. (C-F) Distribution of sampling in the (φ, ψ) space by different one-dimensional importance-sampling calculations using dihedral φ as the transition coordinate. The color of each point reflects the corresponding sample at a certain time. The unit of σ is 2 degrees, consistent with the bin width in this example. To further compare the performance of meta-eABF with that of MtD and eABF, a realistic application, namely the shuttling of a rotaxane in DMSO,37 was studied using the aforementioned algorithms (Figure 3). Here, orthogonal-space-sampling is crucial because the translational movement of the macrocycle is strongly coupled with the desolvation of the molecular wire. A 200-ns simulation without stratification, i.e., using one very large window spanning the entire reaction pathway, was performed for each algorithm. The reference was obtained by an ABF simulation in excess of 1 μs.37 In this reference simulation, the transition coordinate was stratified into 28 small, 1-Å 12


Page 12 of 27

Page 13 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


wide windows, and the force was collected in bins 0.1 Å wide. As shown in Figure 3A, the PMF calculated from a meta-eABF simulation is very close to the reference one. The minute difference between 32 and 35 Å is probably an artefact arising from the use of small windows in the reference ABF simulation. In this particular example, convergence of the meta-eABF simulation required a computer time roughly one fifth of that of the reference simulation. The general shape of the free-energy profile obtained with eABF is similar to that of the reference simulation, but the calculation is conspicuously still far from being converged. In the case of MtD, the calculated PMFs significantly differ from the reference one, no matter what height and width of Gaussians we chose. We suggest that the size of the window is too large for an MtD simulation, and that if large Gaussian heights and widths are used, the molecular assembly is likely to be driven far from thermodynamic equilibrium, leading to incorrect results.

13



Figure 3. Simulations characterizing the shuttling of a rotaxane in DMSO. Upper: the structure of the rotaxane and the reaction coordinate model. The reaction coordinate model is the distance separating the centroid of the cyclodextrin from the left linker. (A) Corresponding one-dimensional PMF calculated by different algorithms. W=0.1 kcal/mol and 2σ=5 bin widths are used in the meta-eABF simulation. (B) Time evolution of the root-mean-square deviation (RMSD) over the free-energy landscape of the reference. To examine the convergence rate of each algorithm, the time evolution of the RMSD of the free-energy profile with respect to the reference one was monitored for each simulation. As depicted in Figure 3B, meta-eABF possesses convergence properties vastly superior to those of both eABF and MtD. Interestingly enough, although plain eABF indicates a moderate convergence rate initially, improvement occurs slowly after 100 ns, because correcting a misevaluated average force with a large weight is admittedly difficult.

The results depicted in Figure 3 underscore a remarkable convergence rate for meta-eABF. Moreover, it should be noted that in ABF-based algorithms, the convergence rate will quadruple if the window range is reduced by half.38 However, meta-eABF obviates the need for a stratification strategy through the introduction of an MtD-like memory kernel, which quickly makes the sampling distribution uniform along the transition coordinate, even in the absence of a staging strategy. The significant efficiency of meta-eABF depicted in Figure 3 clearly shows that this algorithm exploits the best features of MtD and eABF, without their shortcomings and limitations, as mentioned above.

In order to show the performance of meta-eABF for the exploration of

14


Page 14 of 27

Page 15 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


multidimensional free-energy surfaces with a complex, rugged orthogonal space, translocation of a halide ion through a transmembrane synthetic peptide nanotube was examined. Unlike that of small monovalent cations like sodium, transport of chloride by the tubular structure is markedly unfavorable. The halide ion was found in a previous computational investigation to graze the wall periodically during its translocation to compensate for suboptimal hydration.34 Here, a two-dimensional free-energy calculation characterizing the translocation of a chloride ion across a tubular structure consisting of eight stacked cyclo[LW]4 cyclic peptides (where L and W denote, respectively, D-leucine and L-tryptophan) was performed. The distance separating the ion from the center of mass of the nanotube projected along its longitudinal axis, and that in the radial direction were chosen as the reaction-coordinate model, as described in Figure 4A.

Figure 4. Simulations characterizing the translocation of a halide ion through a transmembrane synthetic peptide nanotube. (A) Reaction-coordinate model. Lipid bilayer and water are not shown for clarity. Two-dimensional free-energy landscape obtained by (B) a 20-ns meta-eABF, (C) a 20-ns eABF, and (D) a 200-ns eABF simulation without stratification strategy. W=0.1 kcal/mol and 2σ=5 bin width are 15



used in the meta-eABF simulations. All other parameters of the three simulations are identical. The result of a 20-ns meta-eABF simulation without any stratification strategy is depicted in Figure 4B. Convergence is achieved within roughly one sixth of the time of an ABF simulation.34 It is worth noting in the latter ABF simulation34 the transition pathway was stratified, thereby trading off rigorousness of orthogonal space sampling (Figure S3) for efficiency.38 In the absence of a stratification strategy, the sampling efficiency will be greatly reduced (see the results of a 20-ns eABF simulation in Figure 4C and two 20-ns MtD simulations in Figure S8), indicating that meta-eABF can balance out the efficiency of transition-coordinate-space sampling and reliability of orthogonal-space sampling. For eABF simulation, a significant improvement can be observed by extending the simulation timescale from 20 to 200 ns (see Figure 4D), nevertheless, the height of many peaks is still suboptimally evaluated on the generated free-energy surface. This result clearly demonstrates a superior performance of meta-eABF in mapping rugged multidimensional free-energy surfaces compared with conventional eABF.

The zigzagging lowest free-energy path reflects the diffusion mode of the chloride ion in the hollow tubular structure. To further rationalize grazing of its wall by the monovalent anion, an analysis of the trajectory was performed. Owing to the exiguous hollow cavity, the ion coordination is 4 in the structure depicted in Figure 4A. The coordination number of the anion inside the nanotube is, therefore, much lower than that in bulk water. Since the halide ion has to traverse the nanotube partially naked, it

16


Page 16 of 27

Page 17 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


forms whenever possible transient interactions with the N—H groups of the nanotube to compensate for its incomplete hydration (Figure S7).

Figure 5. Free-energy landscapes describing the translocation of a halide ion through a transmembrane peptide nanotube by means of a 20-ns (A) MtD + ABF estimator (W=0.1 kcal/mol, 2σ=5 bin width), (B) MtD + ABF estimator (W=0.25 kcal/mol, 2σ=14 bin width), (C) eMtD + CZAR estimator (W=0.1 kcal/mol, 2σ=5 bin width) and (D) eMtD + CZAR estimator (W=0.25 kcal/mol, 2σ=14 bin width) simulations without stratification strategy. All other parameters of the simulations are identical. To further clarify the efficiency of meta-eABF, using the same computational assay featuring the peptide nanotube formed by eight stacked cyclo[LW]4 cyclic peptides, we thoroughly compared the performance of meta-eABF with that of all the variants of ABF and MtD available in the Colvars module of NAMD, namely ABF, eABF, MtD,

MtD

plus

an

ABF-like

estimator

(MtD-ABF

hybrid

method),

extended-Lagrangian-based metadynamics (eMtD), eMtD plus an eABF estimator (eMtD-eABF hybrid method), and well-tempered MtD (WTMtD).

As shown in Figure S9 and Figure S10, none of the simulations performed with the traditional importance-sampling algorithms did converge within 100 ns, except for 17



that using WTMtD (W=0.1 kcal/mol, 2σ=5 bin width, T = 4200 K), thereby suggesting that meta-eABF can be as much as five times faster than classical ABF, eABF, MtD and eMtD. WTMtD (W=0.1 kcal/mol, 2σ=5 bin width, T = 4200 K) is admittedly able to recover the correct PMF within 58 ns, but is still about 2.9 times slower than meta-eABF. Introducing the ABF/eABF estimators into MtD/eMtD can both improve the convergence rate and reduce the dependence on the simulation parameters, as depicted in Figure 5. Interestingly enough, those simulations relying upon an MtD + ABF estimator (W=0.25 kcal/mol, 2σ=14 bin width) converged in nearly 20 ns, that is within the same ballpark as meta-eABF. Still, meta-eABF remains faster than any of the latter methods, as can be concluded from Table 1. Moreover, one should keep in mind that use of classical TI and ABF (estimators) requires knowledge of the underlying internal forces. This prerequisite may greatly limit the scope of applications of such gradient-based methods. In stark contrast, meta-eABF, conventional eABF and (e)MtD do not suffer from this limitation.

To highlight the ability of meta-eABF to handle challenging problems, an additional simulation characterizing the permeation of DDA (2’,3’-dideoxyadenosine) across a 2:1 POPC-cholesterol membrane was carried out (Figure S11). The PMF for the translocation of the permeant along the normal to the lipid bilayer, i.e., the z-direction of Cartesian space, obtained from a 510-ns meta-eABF simulation, splitting the reaction pathway into three 15-Å wide, non-overlapping windows, is shown in Figure S11B. The reference free-energy profile was inferred from a 2.467-µs ABF simulation, using nine windows.43 It is apparent that the PMFs obtained with 18


Page 18 of 27

Page 19 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


meta-eABF and conventional ABF are nearly identical, though convergence of ABF is more susceptible to the slow collective movements of the lipids, resulting in longer simulation times. The present series of free-energy calculations indicates that meta-eABF sweeps the entire range of z-values very rapidly, yet without perceptible loss of accuracy.

From the above performance tests in the toy models and practical applications, as summarized in Table 1 and Table S1, it can be concluded that meta-eABF is adaptive, robust and possesses a self-learning ability. We, therefore, suggest that meta-eABF be used as a first choice for any general application.

Table 1. Convergence rate (with respect to meta-eABF) for different importance-sampling algorithms obtained from simulations describing the translocation of a halide ion through a transmembrane peptide nanotube. Methods

Convergence rate

Meta-eABF

1.0

ABF

< 0.2

eABF

< 0.2

MtD (W=0.1 kcal/mol, 2σ=5 bin width) MtD (W=0.25 kcal/mol, 2σ=14 bin width) eMtD (W=0.1 kcal/mol, 2σ=5 bin width) eMtD (W=0.25 kcal/mol, 2σ=14 bin width)

Methods

< 0.2 < 0.2 < 0.2 < 0.2

WTMtD (W=0.1 kcal/mol, 2σ=5 bin width, bias temperature T = 4200 K) WTMtD (W=0.25 kcal/mol, 2σ=14 bin width, T = 4200 K) WTMtD (W=0.1 kcal/mol, 2σ=5 bin width, T = 5000 K) MtD (W=0.1 kcal/mol, 2σ=5 bin width) + ABF estimator MtD (W=0.25 kcal/mol, 2σ=14 bin width) + ABF estimator eMtD (W=0.1 kcal/mol, 2σ=5 bin width) + CZAR estimator eMtD (W=0.25 kcal/mol, 2σ=14 bin width) + CZAR estimator

Convergence rate 0.34 < 0.2 < 0.2 0.43 ~1.0 < 0.25 0.4

Although we did not focus on the efficiency of orthogonal-space sampling in the development of the meta-eABF algorithm, the results of the present work underscore 19



Page 20 of 27

that it represents a significant improvement over both eABF and MtD. Further enhancement of sampling efficiency can be expected, as it is possible to combine seamlessly ABF-based schemes with ergodic-sampling approaches, either through a multiple-walker strategy,39 or using non-equilibrium switches.40 Recently, neural networks have also been associated to importance-sampling algorithms,41,42 which also constitutes an appealing alternative for orthogonal-space sampling.

ASSOCIATED CONTENT

Supporting Information

The Supporting Information is available free of charge on the ACS Publications website at DOI: xxxxx/xxxxxx

Time evolution of ΔA’eABF and ΔA’MtD. A Colvars configuration file showing the usage of meta-eABF in NAMD. The difficulty of orthogonal sampling and the unreliability of a stratification strategy. Simulations describing a particle moving on a model potential energy surface. The effect of hillWeight and hillWidth on meta-eABF simulation. Comparison between meta-eABF and the MtD-ABF hybrid method. Free-energy profile characterizing the shuttling of a rotaxane in DMSO using the

MtD-ABF

hybrid

method.

Comparison

of

meta-eABF

and

other

importance-sampling algorithms. Simulations for the translocation of a halide ion through a transmembrane synthetic peptide nanotube. Discussion about poor choice of

20


Page 21 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


free-energy calculation method. Methodological details. (PDF)

Movies illustrating that poor choice of free-energy calculation method can lead to unreliable results. (ZIP)

AUTHOR INFORMATION Corresponding Author *Email: [email protected] (W. C.), [email protected] (C. C.) Notes The authors declare no competing financial interests. ACKNOWLEDGMENT This work was supported by the National Natural Science Foundation of China (21773125), the Natural Science Foundation of Tianjin, China (18JCYBJC20500), and the Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (Second Phase) under Grant U1501501. The GENCI and CINES, Montpellier, France are gratefully acknowledged for provision of generous amounts of CPU time. C. C. is indebted to the Centre National de la Recherche Scientifique for the support of his joint research program (PICS) with the People’s Republic of China. The authors express their gratitude to Tony Lelièvre (CERMICS, École des Ponts ParisTech) and Benoît Roux (University of Chicago) for insightful discussions.

21



REFERENCES (1) Maximova, T.; Moffatt, R.; Ma, B.; Nussinov, R.; Shehu, A. Principles and overview of sampling methods for modeling macromolecular structure and dynamics. PLoS Comput. Biol. 2016, 12 (4), e1004619. (2) Chipot, C. Frontiers in free-energy calculations of biological systems. WIRES Comput. Mol. Sci. 2014, 4 (1), 71-89. (3) Bernardi, R. C.; Melo, M. C.; Schulten, K. Enhanced sampling techniques in molecular dynamics simulations of biological systems. Biochim. Biophys. Acta 2015, 1850 (5), 872-877. (4) Okamoto, Y. Generalized-ensemble algorithms: enhanced sampling techniques for Monte Carlo and molecular dynamics simulations. J. Mol. Graph. Model. 2004, 22 (5), 425-439. (5) Lelièvre, T.; Stoltz, G.; Rousset, M. Free energy computations: A mathematical perspective. World Scientific: 2010. (6) Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23 (2), 187-199. (7) Kästner, J. Umbrella sampling. WIRES Comput. Mol. Sci. 2011, 1 (6), 932-942. (8) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99 (20), 12562-12566.

22


Page 22 of 27

Page 23 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(9) Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. WIRES Comput. Mol. Sci. 2011, 1 (5), 826-843. (10) Huber, T.; Torda, A. E.; van Gunsteren, W. F. Local elevation: a method for improving the searching properties of molecular dynamics simulation. J. Comput. Aided Mol. Des. 1994, 8 (6), 695-708. (11) Grubmüller, H. Predicting slow structural transitions in macromolecular systems: Conformational flooding. Phys. Rev. E 1995, 52 (3), 2893. (12) Darve, E.; Pohorille, A. Calculating free energies using average force. J. Chem. Phys. 2001, 115 (20), 9169-9183. (13) Comer, J.; Gumbart, J. C.; Hénin, J.; Lelièvre, T.; Pohorille, A.; Chipot, C. The adaptive biasing force method: Everything you always wanted to know but were afraid to ask. J. Phys. Chem. B 2014, 119 (3), 1129-1151. (14) Darve, E.; Rodríguez-Gómez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations. J. Chem. Phys. 2008, 128 (14), 144120. (15) Zheng, L.; Yang, W. Practically efficient and robust free energy calculations: double-integration orthogonal space tempering. J. Chem. Theory Comput. 2012, 8 (3), 810-823. (16) 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.

23



2008, 105 (51), 20227-20232. (17) Fu, H.; Shao, X.; Chipot, C.; Cai, W. Extended adaptive biasing force algorithm. An on-the-fly implementation for accurate free-energy calculations. J. Chem. Theory Comput. 2016, 12 (8), 3506-3513. (18) Lesage, A.; Lelièvre, T.; Stoltz, G.; Hénin, J. Smoothed biasing forces yield unbiased free energies with the extended-system adaptive biasing force method. J. Phys. Chem. B 2017, 121 (15), 3676-3685. (19) Zhao, T.; Fu, H.; Lelièvre, T.; Shao, X.; Chipot, C.; Cai, W. The extended generalized adaptive biasing force algorithm for multidimensional free-energy calculations. J. Chem. Theory Comput. 2017, 13 (4), 1566-1576. (20) Comer, J.; Phillips, J. C.; Schulten, K.; Chipot, C. Multiple-replica strategies for free-energy calculations in NAMD: multiple-walker adaptive biasing force and walker selection rules. J. Chem. Theory Comput. 2014, 10 (12), 5276-5285. (21) Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100 (2), 020603. (22) Tiwary, P.; Parrinello, M. A Time-independent free energy estimator for metadynamics. J. Phys. Chem. B 2015, 119 (3), 736-742. (23) Piana, S.; Laio, A. A bias-exchange approach to protein folding. J. Phys. Chem. B 2007, 111 (17), 4553-4559.

24


Page 24 of 27

Page 25 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(24) Marinelli, F.; Pietrucci, F.; Laio, A.; Piana, S. A kinetic model of trp-cage folding from multiple biased molecular dynamics simulations. PLoS Comput. Biol. 2009, 5 (8), e1000452. (25) Chen, M.; Cuendet, M. A.; Tuckerman, M. E. Heating and flooding: A unified approach for rapid generation of free energy surfaces. J. Chem. Phys. 2012, 137 (2), 024102. (26) Cuendet, M. A.; Tuckerman, M. E. Free energy reconstruction from metadynamics or adiabatic free energy dynamics simulations. J. Chem. Theory Comput. 2014, 10 (8), 2975-2986. (27) Babin, V.; Roland, C.; Sagui, C. Adaptively biased molecular dynamics for free energy calculations. J. Chem. Phys. 2008, 128 (13), 134101. (28) Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient exploration of reactive potential energy surfaces using Car-Parrinello molecular dynamics. Phys. Rev. Lett. 2003, 90 (23), 238302. (29) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26 (16), 1781-1802. (30) Fiorin, G.; Klein, M. L.; Hénin, J. Using collective variables to drive molecular dynamics simulations. Mol. Phys. 2013, 111 (22-23), 3345-3362.

25



(31) Chen, H.; Fu, H.; Shao, X.; Chipot, C.; Cai, W. ELF: An extended-lagrangian free energy calculation module for multiple molecular dynamics engines. J. Chem. Inf. Model. 2018, 58(7), 1315-1318. (32) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185 (2), 604-613. (33) Mones, L.; Bernstein, N.; Csányi, G. Exploration, sampling, and reconstruction of free energy surfaces with gaussian process regression. J. Chem. Theory Comput. 2016, 12 (10), 5100-5110. (34) Hé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 (1), 35-47. (35) Wang, Z. X.; Duan, Y. Solvation effects on alanine dipeptide: A MP2/cc-pVTZ//MP2/6-31G** study of (Φ, Ψ) energy maps and conformers in the gas phase, ether, and water. J. Comput. Chem. 2004, 25 (14), 1699-1716. (36) Chipot, C.; Pohorille, A. Conformational equilibria of terminally blocked single amino acids at the water−hexane interface. A molecular dynamics study. J. Phys. Chem. B 1998, 102 (1), 281-290. (37) Liu, P.; Chipot, C.; Shao, X.; Cai, W. Solvent-controlled shuttling in a molecular switch. J. Phys. Chem. C 2012, 116 (7), 4471-4476.

26


Page 26 of 27

Page 27 of 27 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(38) Comer, J.; Roux, B.; Chipot, C. Achieving ergodic sampling using replica-exchange free-energy calculations. Mol. Simul. 2014, 40 (1-3), 218-228. (39) Minoukadeh, K.; Chipot, C.; Lelievre, T. Potential of mean force calculations: a multiple-walker adaptive biasing force approach. J. Chem. Theory Comput. 2010, 6 (4), 1008-1017. (40) Suh, D.; Radak, B. K.; Chipot, C.; Roux, B. Enhanced configurational sampling with hybrid non-equilibrium molecular dynamics–Monte Carlo propagator. J. Chem. Phys. 2018, 148 (1), 014101. (41) Guo, A. Z.; Sevgen, E.; Sidky, H.; Whitmer, J. K.; Hubbell, J. A.; Pablo, J. J. d. Adaptive enhanced sampling by force-biasing using neural networks. J. Chem. Phys. 2018, 148 (13), 134108. (42) Galvelis, R.; Sugita, Y. Neural network and nearest neighbor algorithms for enhancing sampling of molecular dynamics. J. Chem. Theory Comput. 2017, 13 (6), 2489-2500. (43) Tse, C. H.; Comer, J.; Wang, Y.; Chipot, C. The link between membrane composition and permeability to drugs. J. Chem. Theory Comput. 2018, 14 (6), 2895−2909.

27


Zooming Across the Free-Energy Landscape, Shaving Barriers and

Recommend Documents