The Extended Generalized Adaptive Biasing Force Algorithm for

Subscriber access provided by UB + Fachbibliothek Chemie | (FU-Bibliothekssystem)

Article

The Extended Generalized Adaptive Biasing Force Algorithm for Multidimensional Free-Energy Calculations Tanfeng Zhao, Haohao Fu, Tony Lelièvre, Xueguang Shao, Christophe Chipot, and Wensheng Cai J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00032 • Publication Date (Web): 02 Mar 2017 Downloaded from http://pubs.acs.org on March 6, 2017

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 free 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 accessible to all readers and 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.

Journal of Chemical Theory and Computation 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 36

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

Journal of Chemical Theory and Computation

The Extended Generalized Adaptive Biasing Force Algorithm for Multidimensional Free-Energy Calculations Tanfeng Zhao,†,┴ Haohao Fu,†,┴ Tony Lelièvre,∇ 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 (Nankai University), Tianjin 300071, China

∇Université

$

Paris-Est, CERMICS (ENPC), INRIA, 77455 Marne-la-Vallée, France

Laboratoire International Associé Centre National de la Recherche Scientifique et University of

Illinois at Urbana-Champaign, Unité Mixte de Recherche No. 7565, Université de Lorraine, B.P. 70239, 54506 Vandœuvre-lès-Nancy cedex, France §

Theoretical and Computational Biophysics Group, Beckman Institute, University of Illinois at

Urbana-Champaign, Urbana, Illinois 61801, United States ¶

Department of Physics, University of Illinois at Urbana−Champaign, 1110 West Green Street,

Urbana, Illinois 61801, United States

* To whom correspondence should be addressed ┴ Contributed equally to this work

1

ACS Paragon Plus Environment


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 36

ABSTRACT Free-energy calculations in multiple dimensions constitute a challenging problem, owing to the significant computational cost incurred to achieve ergodic sampling. The generalized adaptive biasing force (gABF) algorithm calculates n one-dimensional lists of biasing forces to approximate the n-dimensional matrix by ignoring the coupling terms ordinarily taken into account in classical ABF simulations, thereby, greatly accelerating sampling in the multidimensional space. This approximation may however occasionally lead to poor, incomplete exploration of the conformational space compared to classical ABF, especially when the selected coarse variables are strongly coupled. It has been found that introducing extended potentials coupled to the coarse variables of interest can virtually eliminate this shortcoming, and, thus, improve the efficiency of gABF simulations. In the present contribution, we propose a new free-energy method, coined extended generalized ABF (egABF), combining gABF with an extended Lagrangian strategy. The results for three illustrative examples indicate that (i) egABF can explore the transition coordinate much more efficiently compared with classical ABF, eABF and gABF, in both simple and complex cases, and (ii) egABF can achieve a higher accuracy than gABF, with a root mean-squared deviation between egABF and eABF free-energy profiles on the order of kBT. Furthermore, the new egABF algorithm outruns the previous ABF-based algorithms in high-dimensional free-energy calculations, and, hence, represents a powerful importance-sampling alternative for the investigation of complex chemical and biological processes.

2


Page 3 of 36

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


INTRODUCTION Brute-force molecular dynamics (MD) simulations regarded as a “computational microscope” have been widely used to explore biological processes at the atomic level.1-2 In some cases, however, due to the inherent complexity of the biological objects at hand, the computational effort incurred in covering the relevant timescales is incompatible with classical MD simulations run on common computer architectures. Importance sampling, a set of methods that encourage sampling along arbitrarily chosen coarse variables, is one of the most popular avenues to accelerate exploration of the participating slow degrees of freedom.3 Although of a different nature, adaptive biasing force (ABF),4-6 metadynamics,7-9 and umbrella sampling10-11 are commonly used importance-sampling methods . In practice, however, many processes cannot be properly characterized by one, naive coarse variable based on chemical intuition.12-15 Multidimensional free-energy calculations along a set of coarse variables should, therefore, be applied in these cases. To improve ergodic sampling and, thus, reduce the significant computational cost of such simulations, many strategies aimed at improving the convergence of importance-sampling algorithms have been developed. For example, a multiple-walker strategy was combined with umbrella sampling,16-17 metadynamics18 and ABF19-20 to improve the convergence rate of the importance-sampling scheme. The self-learning adaptive umbrella sampling21-22 and the on-the-path random walk methods23 focus on sampling along the least free-energy pathway, which is shown to diminish the computational cost of the simulation. Although it is not an importance-sampling algorithm, accelerated MD,24-25 embodies the idea of trading accuracy for efficiency in free-energy calculations.26 To enhance ergodic sampling of ABF-based methods, Chipot and Lelièvre have proposed the generalized ABF (gABF) algorithm.27 This method can greatly accelerate sampling in free-energy

3



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 4 of 36

calculations in multidimensional space by ignoring the coupling terms inherent to classical ABF.27 A qualitatively correct multidimensional free-energy landscape can be obtained within a fraction of the time required to reach convergence of the time-dependent biases in ABF simulations. The approximate free-energy landscape supplied by gABF can be used as an initial guess to expedite a classical ABF calculation, from whence accurate free energies are determined. Two limitations of the above strategy, however, greatly thwart its widespread application. First, the chosen coarse variables must obey the prerequisites of classical ABF.28 Second, on account of the omission of the coupling terms, the bias obtained from a gABF calculation under certain circumstances deviates markedly from the free-energy landscape, thereby limiting the efficiency of the algorithm to explore the available conformational space.27 It was found, however, that the extended ABF (eABF) method3, 28-31 effectively decouples the coarse variables utilized to model the transition coordinate, by introducing extended degrees of freedom. Due to the orthogonality of the extended variables, the contribution of the corresponding coupling terms is reduced considerably. Taking into account both accuracy and efficiency, the main thrust of the present work is to incorporate an extended-Lagrangian strategy into gABF. By combining gABF with eABF,28 we propose a novel approach, coined extended generalized ABF (egABF), expected to be able of eliminating the requirements of ABF and improve both the accuracy and efficiency of geometrical free-energy calculations. We have developed a code to perform egABF calculations based on the Colvars module,32 obviating the need of a correction of the gABF simulation followed by standard ABF.

4


Page 5 of 36

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


METHODS Theoretical background Detailed theoretical background about ABF, eABF, gABF and egABF can be found in the Appendix A. Here, we briefly outline the idea of the new egABF algorithm. The overdamped Langevin dynamics of standard ABF27 can be expressed as,

(1)

where function,

denotes an n-dimensional transition coordinate,

is the potential energy

is the so-called instantaneous collective force (or local mean force) associated with

(see also equation (8) in the Appendix B),

is the biasing function associated with

Wiener process or white noise that underlies the random force, and

, with

,

is the being the

temperature. Originally, gABF was proposed to address the overwhelming computational cost of free-energy calculations relying on an n-dimensional transition coordinate, with n ≥ 2. The main idea of gABF consists in replacing the n-dimensional bias function, n one-dimensional functions

, with the sum of

(see equation 12 in the Appendix B). In other words, the coupling

terms between each coarse variable are neglected. Following this strategy, it is, therefore, apparent that mapping the free-energy landscape would require far less computational time. Under the assumption that the coarse variables, i.e.,

, are completely decoupled, gABF will generate a biasing

force nearly identical to that of standard ABF. To recover the free-energy landscape, one needs to compute

by means of, ,

5


(2)


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

where

Page 6 of 36

is composed of n biasing forces associated to all the coarse variables.

In practice, the premise that the transition coordinate consists of decoupled coarse variables is by and large never verified, which actually greatly limits the direct application of the gABF algorithm. As a possible solution to this problem, the eABF algorithm helps decouple the coarse variables by introducing an extended potential to the system of interest.28 In eABF free-energy calculations, a biasing force is applied to a fictitious particle, which is coupled to the coarse variable by means of a stiff harmonic spring, , where

is the force constant of the spring,

(3)

is a fictitious degree of freedom.

In this work, we propose to use the gABF approach on the extended system of the eABF algorithm, which yields the egABF dynamics,

(4)

where

is an n-dimensional Wiener process,

is a one-dimensional function of

is the coarse variable, and

. Our primary objective is to combine the advantages of the eABF

and gABF algorithms outlined above. Choosing a sufficiently small value of

for the extended

degrees of freedom, the coupling of the connected coarse variables will become marginal (see Appendix B for detailed discussion). Under this condition, the assumption of small coupling for an accurate gABF calculation is satisfied. In other words, the sum of n one-dimensional biasing functions of the gABF algorithm approximates to the n-dimensional biasing function in the ABF algorithm. However, an efficient exploration of the transition coordinate in the eABF algorithm requires a

6


Page 7 of 36

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


relatively large value of

to guarantee the accuracy of the calculated free-energy change.28 Proper

choice of the force constant is, therefore, important when using the egABF algorithm, as will be seen in the Result and Discussion section, as well as in the Appendix B.

Implementation Details To carry out a gABF/egABF calculation using the popular molecular dynamics engine NAMD,33 we have implemented an on-the-fly code interfaced with the Colvars32 module through its Colvars scripting feature. Our implementation is readily applicable to a hybrid QM/MM (quantum mechanics/molecular mechanics) simulation to investigate complex chemical reaction mechanisms. In egABF calculations, the on-the-fly eABF framework developed in our previous work28 is called at every time step. This implementation requires a minimum amount of human interventions. Only a number of parameters in the NAMD configuration file are needed to specify the frequency to output data and the names of the input and output files (see the Supporting Information).

Free-Energy Calculations As shown in Figure 1, three molecular assemblies were selected to demonstrate the applicability, the efficiency, and the accuracy of the new method. For comparison purposes, classical ABF, eABF, gABF, and egABF were employed to determine the two-dimensional free-energy landscape characterizing the isomerization of N-acetyl-N’-methylalanylamide (NANMA), which has been recurrently utilized to benchmark the performance of new free-energy methods.19,20,32,34-38 Next, eABF and egABF were used to carry out the free-energy calculations describing the complex membrane translocation of aspirin39,40 and the isomerization of deoxycytidine.41,42 The coarse variables depicted

7



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 8 of 36

in Figure 1C are highly coupled, as shown in our previous work.28 The least free-energy pathways (LFEP) connecting the minima representing the pseudo-rotation of furanose in deoxycytidine on two-dimensional free-energy landscapes were identified using the LFEP algorithm.43 Based on the guidelines that we put forth for the choice of the parameters in eABF calculations,28 the parameter used to determine

of the spring, i.e., the standard deviation between the fictitious particle and the

associated coarse variable, was equal to the bin width of the collective variable, and the oscillation period of the fictitious particle was set to 200 fs.

Figure 1. Molecular assemblies and the corresponding transition coordinates used in this study. (A) Isomerization of N-acetyl-N’-methylalanylamide (NANMA) in water. The transition coordinate consists of the two backbone dihedral angles,

and

transition coordinate is defined as a torsional angle

. (B) Membrane translocation of aspirin. The and the distance d separating the center of mass

of the solute from that of the lipid bilayer (determined from its heavy atoms) projected onto z axis. (C) Isomerization of deoxycytidine. The transition coordinate is defined as the linear combination of the two dihedral angles,

and

.

8


Page 9 of 36

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


Simulation Details In this study, all the MD simulations were carried out using the parallel, scalable program NAMD 2.1133 with the CHARMM force field,44,45 and the CHARMM general force field (CGenFF)46,47. The VMD program was used to display and analyze the MD trajectories.48 The deoxycytidine was placed in vacuum. A time step of 0.5 fs was utilized to integrate the equations of motion. The temperature was maintained at 300 K employing Langevin dynamics.49 The NANMA and the membrane system with an aspirin molecule were solvated in periodic boxes of TIP3P50 water molecules. For them, the r-RESPA multiple time-step algorithm51 was employed to integrate the equations of motion with a time step of 1 and 2 fs for short- and long-range interactions, respectively. Covalent bonds involving hydrogen atoms of (bio)organic molecules were constrained to their equilibrium length by means of the SHAKE/RATTLE52,53 algorithms while the SETTLE54 algorithm was used for water molecules. The temperature and the pressure were maintained at 300 K and 1 atm, respectively, using Langevin dynamics and the Langevin piston method.49 Long range electrostatic forces were evaluated by means of the particle mesh Ewald algorithm.55 The short-range van der Waals and electrostatic interactions were calculated by the smoothed 12.0 Å spherical cutoff.

ab initial Calculation For comparison purposes, the absolute free energy of the stable C2’-endo and C3’-endo states of deoxycytidine was determined quantum-mechanically. The free-energy difference between these two states was used to evaluate the accuracy of the egABF calculation. The initial conformations were taken from the molecular dynamics trajectory. The G4(MP2)56 composite approach was used. All the quantum-mechanical calculations were performed, employing the Gaussian 09 package.57

9



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 10 of 36

RESULTS AND DISCUSSION Accuracy and Efficiency of egABF In this contribution, the isomerization process of NANMA was investigated by means of the classical ABF, eABF, gABF and egABF algorithms in order to assess the performance of the new egABF scheme. The free-energy landscape reproduced from a 200 ns ABF calculation was used as a reference, as shown in Figure S1 of the Supporting Information. Infrared and Raman spectra experiments37,38 show that the αR and β conformations (labeled in Figure 3B) are the two major conformations observed in aqueous solution, congruent with our results highlighted in Figure 3B.

10


Page 11 of 36

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


Figure 2. Time evolution of the two-dimensional free-energy landscape underlying the isomerization of NANMA, obtained from classical ABF (A-D), gABF (E-H), eABF (I-L), and egABF (M-P) calculations. The RMSD of each PMF with respect to a reference one derived from a 200 ns ABF calculation is mentioned under each panel.

Figure 2 shows the free-energy maps charactering the isomerization of NANMA, obtained by means of the four methods with progressively increasing simulation time. The root-mean-square

11



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 12 of 36

deviation (RMSD) of the free-energy landscape with respect to the ABF reference is calculated at different stages of the simulations (see Figure S1 in the Supporting Information). By comparing the time evolution of the different RMSDs in Figure 2, it can be seen that (i) gABF is able to locate the minima well within a very short simulation time, (ii) although eABF does not reproduce the free-energy landscape very well in the early stage of the simulation, it rapidly refines it subsequently, (iii) by combining the advantages of both gABF and eABF, egABF always yields the smallest RMSD at each stage, compared to the other three methods. The free-energy landscapes generated from two long gABF and egABF calculations are depicted in Figure 3. The RMSD of the free-energy landscape determined from the egABF simulation with respect to the ABF reference is decreased to 0.6 kcal/mol, which is much lower than that corresponding to the gABF simulation. The basins and peaks in Figure 3B resemble the ones observed in classical ABF simulations (Figure S1), while those in Figure 3A are still slightly distorted compared to the reference free-energy landscape. The difference between the free-energy landscapes of the converged gABF and ABF simulations clearly illuminates the loss of information caused by ignoring the coupling terms accounted for in ABF simulations. The eABF method, however, can alleviate this shortcoming by means of its extended potential, which decouples the coarse variables. Put together, egABF appears to be more efficient than both eABF and classical ABF, and, in practice, more accurate than gABF.

12


Page 13 of 36

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


Figure 3. Two-dimensional free-energy landscapes underlying the isomerization of NANMA, obtained from two 120 ns (A) gABF, and (B) egABF calculations. RMSD of each PMF with respect to that by a 200 ns ABF calculation is listed under each panel.

Membrane Translocation of Aspirin The second application of egABF is the complex passive diffusion of aspirin across a membrane. The two-dimensional free-energy landscapes underlying this process obtained from different time of eABF and egABF calculations, are depicted in Figure 4. A previous study40 suggests that the hydrogen of the carboxyl group exhibits a strong propensity towards a cis conformation in the membrane, because of the 1,4 intramolecular interaction of the hydroxyl hydrogen and the carbonyl oxygen. In the glycerol region of membrane, other than the 1,4 interaction, the carboxyl group of aspirin in trans conformation can form hydrogen bond with the esters of the lipids, which means that aspirin has two stable states at the water/membrane interface.

13



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 14 of 36

Figure 4. Time evolution of two-dimensional free-energy landscape underlying the membrane translocation of aspirin, obtained from (A) eABF, and (B) egABF calculations. Water: d > 20 Å, water/membrane interface: d = 10 - 20 Å, membrane midplane: d = 0 Å.

As may be seen in Figure 4A-D, in a 160–ns eABF simulation, no clear local minimum could be found in the region corresponding to the trans conformation of the carboxyl group (θ = ±180°), due to insufficient sampling. Conversely, the two-dimensional free-energy landscape obtained from a 20-ns egABF simulation features all the basins and peaks emerging at the correct locations, as depicted in Figure 4E. After 160 ns, all the regions of conformational space were well sampled, as shown in Figure 4H. The calculated free-energy difference between the bound and the free state of aspirin (-2.9 kcal/mol) is in agreement with the experimental binding free energy (-3.4 kcal/mol).39,40 The unparalleled sampling efficiency of egABF convincingly suggests that this new method is eminently suited to complex transition coordinates.

14


Page 15 of 36

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


Pseudorotation of Furanose in deoxycytidine To explore further the applicability of egABF to two-dimensional free-energy calculations in which the two coarse variables are highly coupled, a simulation describing the isomerization of deoxycytidine was performed, as depicted in Figure 5B. The transition coordinate was defined as the linear combination of two dihedral angles.42 (5)

(6) For comparison purposes, an eABF simulation was also carried out. The PMFs shown in Figure 5A and 5B clearly describe the mode of isomerization of deoxycytidine, i.e. pseudo rotation, with two stable states, C2’-endo and C3’-endo, the former being more stable than the latter, in accordance with ab initio42,58 and experimental data.41 The free-energy difference between the two conformations was calculated to be 0.8 kcal/mol by eABF, and 0.5 kcal/mol by egABF. Both results agree well with the experimental data (0.3 kcal/mol),41 the quantum-mechanical estimate reported in the literature (0.8 kcal/mol)58 and our own estimate (0.7 kcal/mol). Applying the free-energy landscape shown in Figure 5A as a reference, the RMSD between the free-energy landscapes of figures 5B and 5A was estimated to be 0.6 kcal/mol. The transition pathways connecting the two basins, i.e., the pseudorotation of the five-member ring,42 were determined using the LFEP algorithm.39 The free-energy changes along the putative minimum-action paths are shown in Figure 6. The two profiles appear to be very close to each other, suggesting that the PMF calculated by egABF can characterize the pseudorotation of a furan ring with appropriate accuracy. Our findings further suggest that egABF is perfectly suited to those cases wherein the coarse variables are highly coupled.

15



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 16 of 36

Figure 5. The free-energy landscapes charactering the isomerization of deoxycytidine, obtained from (A) 30 ns eABF and (B) 30 ns egABF simulations.

Figure 6. Free-energy profiles as a function of the position (s) along the LFEPs. s is the order parameter of the pathway. s = 0.0 and s = 1.0 correspond to the centers of the basins on the left- and right-hand side of each free-energy map.

16


Page 17 of 36

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


Choice of the Parameters of egABF Calculations As mentioned previously in the Method section, decoupling of the coarse variables in the egABF algorithm is strongly related to the stiffness of the spring between the fictitious particle and the transition coordinate. Although a low force constant is needed to decouple these variables, a relative large stiffness is required in eABF simulations to guarantee a good exploration of the transition coordinate.28 Here, we have explored the effect of the force constant of the spring in egABF simulations. In the Colvars module of NAMD,32 the force constant of the spring,

is calculated as

follows, (7) where σ is the standard deviation between the extended degree of freedom and the associated collective variable. The end-users can specify the value of σ with the keyword extendedFluctuation. To investigate the effect of the spring constant in egABF simulations, the free-energy change characterizing the isomerization of deoxycytidine was determined using different values of σ. The bin width of the collective variables was set in all simulations to 5º. As shown in figures 5B and 7, the free-energy profiles vary drastically with different values of extendedFluctuation, indicating that egABF is acutely sensitive to the stiffness of the spring. The best result compared to the reference two-dimensional free-energy landscape of Figure 5A was obtained when extendedFluctuation was set to the value of the bin width. These results provide a guideline for choosing the parameters of egABF simulations, which coincides with those of eABF.28

17



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 18 of 36

Figure 7. Free-energy landscapes characterizing the isomerization of deoxycytidine, using different values for the egABF parameter, extendedFluctuation.

CONLUSION The strength of adaptive biasing methods compared to other non-adaptive importance-sampling schemes is that the added energy is adapted on the fly in order to flatten the potential along the transition coordinate. ABF-based methods are, therefore, much more convenient to utilize in practice since the end-user does not have to define the whole biasing potential, but only to choose the transition coordinate. Moreover, adaptive biasing force techniques are particularly efficient compared to adaptive

18


Page 19 of 36

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


biasing potential methods, since they directly compute the derivative of the free energy, and not the free energy itself. The derivative is what is actually needed in order to update the Cartesian force in the dynamics. In the event where the coarse variables utilized to model the transition coordinate are highly coupled, the gABF algorithm is able to speed up sampling by ignoring the coupling terms that are traditionally accounted for in classical ABF simulations. The bias computed by the algorithm is, however, far from the actual n-dimensional free energy, so that in all likelihood, the system will remain trapped in some metastable states. The eABF algorithm decouples the coarse variables through the introduction of fictitious particles attached to them by means of stiff springs. The features of gABF and eABF are combined into a novel method coined egABF. Inheriting the high speed of exploration of the transition coordinate model characteristic of gABF, egABF can boost up the convergence rate compared to the classical ABF and eABF algorithms. Even in those cases where the transition coordinate is formed by highly coupled coarse variables, egABF can yield a reasonably precise reproduction of the free-energy landscape over appreciably short simulation times. In the cases of NANMA and deoxycytidine, both the RMSDs of the final free-energy landscapes with respect to their reference are estimated to be on the order of

. These results are suggestive that the newly

implemented egABF algorithm is well suited to the investigation of complex chemical and biological processes. Another noteworthy feature of gABF is its ability to identify the approximate minimum free-energy transition pathway connecting two metastabilities within limited simulation times, as has been demonstrated previously.27 As an avatar of the gABF algorithm, egABF inherits this feature and can markedly improve the accuracy of the PMF following the maximum-likelihood path. Since the key

19



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 20 of 36

features of a free-energy landscape are the basins and the minimum free-energy paths connecting them, egABF appears as a very promising and efficient approach for the exploration of the LFEP in a multidimensional free-energy hyperplane.

ASSOCIATED CONTENT Supporting Information Configuration files for running the gABF/egABF code in NAMD. Two-dimensional free-energy surface characterizing the isomerization of NANMA, obtained from a 200 ns ABF calculation (PDF). egABF code with an example (ZIP). This information is available free of charge via the Internet at http://pubs.acs.org.

AUTHOR INFORMATION Corresponding Authors *

E-mail: [email protected] (W.C.).

*

E-mail: [email protected] (C.C.).

Notes The authors declare no competing financial interest.

Acknowledgements This study is supported by National Natural Science Foundation of China (No. 21373117). The GENCI and CINES, Montpellier, France, and Special Program for Applied Research on Super

20


Page 21 of 36

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


Computation of the NSFC-Guangdong Joint Fund (the second phase) of China are gratefully acknowledged for provision of generous amounts of CPU time. The work of T. Lelièvre is supported by the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) /ERC Grant Agreement number 614492.

Appendix A: Comparison of the ABF, eABF, gABF, and egABF Methods These four methods use the same basic idea of biasing the system with external forces to accelerate sampling of the conformational space. The differences lie in the way the forces are calculated and applied. For clarity, in Scheme 1, the ABF-based algorithms are shown in a flow chart accompanied by the detailed implementation of each method. The accuracy and convergence rate of these methods for the three test cases of this study are illustrated in Scheme 2. Compared to ABF and gABF, eABF and egABF do not have any prerequisite on the transition coordinate, which makes them suitable for a broader range of applications.28

Scheme 1. Flow chart of ABF-based algorithm and the detailed implementation of the four methods.

21



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 22 of 36

Scheme 2. Performance of the four methods in the three test examples in this study.

Appendix B: Mathematical Background 1. The egABF algorithm For the sake of simplicity, we present the egABF algorithm applied to the overdamped Langevin dynamics in dimension d,

though the algorithm readily applies to the phase-space Langevin dynamics. Here, d-dimensional stochastic process,

is a d-dimensional Wiener process, and

is a , with

being the temperature. This dynamics is ergodic with respect to the Boltzmann Gibbs measure,

with

.

Let us first recall the original ABF method.4,

6, 35

For a n-dimensional transition coordinate,

(with n < d), it consists in simulating the dynamics embodied in equation (1). One possible formula for

is:

22


Page 23 of 36

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


(8) where

denotes the (α,γ)-component of the inverse of the Gram matrix, G (which is assumed to

be invertible) defined by, .

(9)

The key idea of the ABF algorithm is that the biasing function,

defined in

equation (1), is an approximation of the gradient of the free energy. This approximation

is built

using the configurations visited up to time . We recall that the free energy is defined by: . In the long-time limit, along

(namely for all

converges to

, so that the force

,

(10) in (1) is zero

). Sampling along the transition

coordinate is, therefore, enhanced. The convergence of the ABF method has been analyzed in various previous contributions,59, 60 and these analyses explain the importance of the choice of the transition coordinate in order to augment the numerical efficiency of the algorithm — the timescale associated with the sampling of the conditional measures

should be much smaller than the

timescale associated with the sampling of the original measure,

.61, 62

egABF is based on the combination of two ideas, which have been proposed to expand the applicability of the ABF algorithm, namely, (i) The extended-ABF idea,3,28 which consists in considering an extended system by adding additional degrees of freedom to the original problem, and the associated extended potential, . Standard ABF is then applied using the simple transition coordinate is available in the Colvars module32 of the NAMD molecular dynamics engine.33

23


(11) . This variant


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 24 of 36

(ii) The generalized-ABF idea,27 which is a variant of ABF, enabling the use of multidimensional transition coordinates. It consists in considering, instead of (1), the following variant,

.

It should be noted that in equation (1), the biasing functions

depend on

(12)

, which makes

them difficult to approximate when n > 3, owing to the large number of degrees of freedom involved in the approximation (which typically scales exponentially with respect to n). In contrast, in equation (12), each function function of

only depends on a one-dimensional parameter (

is a

), which makes the biasing functions easier to approximate.

The common denominator of these two approaches is ABF, an adaptive importance sampling algorithm: One does not try to bias the dynamics by the real free energy associated with , but only by an approximation of it. In eABF, the bias, indeed, converges to an approximation, energy,

, of the free

, defined by a convolution of the associated densities of states, (13)

where ∗ denotes the convolution, and

is the Gaussian kernel with variance

. In gABF, the bias

converges to a tensor product approximation of the free energy, , where

is a primitive of

(14)

. For both eABF and gABF, since the biases can be

shown to converge to a limiting function in the long-time limit (at least under some appropriate assumptions, see below), some standard unbiasing procedures can be used in order to sample the original canonical measure, and, thus, recover the real free energy. The relative efficiency of alternative free-energy estimators27, 28, 63 has been the object of recent discussions in the literature.

24


Page 25 of 36

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 advantage of eABF lies in the simplicity to compute the instantaneous collective forces (they do not involve derivatives of the function , and do not require some non-degeneracy condition on such as det G ≠ 0, see equation (9)). This is particularly useful for situations where holonomic constraints are applied, in which case standard ABF requires to use instantaneous collective forces compatible with these constraints.28 Moreover, a very nice feature of eABF is that, thanks to the convolution procedure of equation (13), the target free energy, energy,

(the smaller

, the smoother

, is smoother than the original free

). In other words,

is typically easier to

approximate than the original target, i.e., larger bins can be used (typically of size less samples are required to get reasonable approximations of chosen too small, in which case

would be far from

. Evidently,

, and the biasing by

) and, thus, should not be would not be

sufficiently efficient to escape from the free-energy wells. In practice, recent numerical experiments28, 64

have shown that in some practical cases of interest, there exists a parameter

that yields an

improved convergence over the original ABF method. The advantage of the gABF approach lies in its ability to deal with large dimensional transition coordinates. Indeed, with standard ABF, the complexity of the method scales exponentially with the number of coarse variables, which, in practice, makes the approach intractable for n > 3. The number of bins becomes so large so that it is impossible to get sufficient statistics in all of them. gABF consists in postulating that the tensor product approximation in equation (14) is sufficiently good to yield a reasonable bias, which helps the system exit free-energy wells. It ought to be noted that a common feature of the two numerical approaches is to reduce the variance of the estimated bias. For eABF, variance reduction is possible thanks to the smoothing induced by the convolution in equation (13). For gABF, variance reduction is a consequence of the fact

25



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 26 of 36

that each sample contributes to the updating of the n one-dimensional functions

. Reducing the

fluctuations on the computed bias helps converge faster to equilibrium.63 In this work, we propose to use gABF on the extended system of the eABF algorithm, yielding the egABF dynamics, which writes,

(15)

where

is an n-dimensional Wiener process,

is a one-dimensional function of

is the coarse variable, and

. The hope is to combine the benefits of the eABF and gABF

approaches outlined above. As discussed in the next section, we expect in particular to improve the long-time convergence of egABF, by comparison with gABF alone, by virtue of a much better decoupling of the coarse variables.

2. Discussion on the longtime convergence In a previous work,27 we analyzed the convergence to equilibrium of gABF. We have shown in a simple two-dimensional example (n = 2, and

) that convergence is exponential with rate

with , where, for

= 1,2,

is the logarithmic Sobolev constant of the family of conditional measures

and .

26


Page 27 of 36

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


Here,

denotes the level sets of

tangential gradient on

, and

is the

,

is the

-th instantaneous collective force defined by equation (8).

It should be stressed that our results hold under an assumption of sufficiently small coupling, namely, . Intuitively, κα measures the coupling of the coordinates orthogonal to

: if

and the other coordinates

, one can check that

(where

) is zero if and only if the dynamics on coordinate dynamics on coordinates the variables

. The parameters

at a fixed value of

is independent from the

measure the mixing property in potential

: the larger

, the faster the components

for

.

Let us now transpose our result on the convergence of gABF to the egABF dynamics described in equation (15). A direct corollary is the following theorem, where we consider for simplicity the following setting — periodic boundary conditions and n = 2. For generalization to other boundary conditions and n > 2, we refer to the original gABF paper.27

Theorem 1 Let us assume that ξ is with values in the 2-dimensional torus in the d-dimensional torus. The potential function

, and that

lives

is assumed to be continuously

differentiable. Then, there exists a stationary measure for the dynamics (15), denoted by . This measure has uniform marginal laws along

Let us denote

, and respectively

and

, the logarithmic Sobolev inequality constants of the

conditional measures with densities, ,

:

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

Page 28 of 36

, respectively. Let us assume that ,

(16)

where . Then, the law of

converges exponentially fast to

in the limit

, with rate

, where . In terms of longtime convergence, the advantage of egABF is that parameter

can be tuned so

that the assumption embodied in equation (16) is satisfied, and the assumption of small coupling, which appears in the original gABF algorithm can be circumvented. In practice,

should be chosen

sufficiently small so that the assumption of equation (16) holds, albeit not too small to prevent the variance

from becoming too large, which, in turn, would deteriorate the approximation of the

free energy through the convolution procedure of equation (13). Our numerical simulations in the main text show that there is, indeed, an appropriate choice for parameter

, which yields a much better

convergence for egABF than for the three other schemes (ABF, eABF and gABF), as illustrated in the various test cases investigated herein.

References (1) Lee, E. H.; Hsin, J.; Sotomayor, M.; Comellas, G.; Schulten, K. Discovery Through the Computational Microscope. Structure 2009, 17, 1295–1306.

28


Page 29 of 36

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


(2) Dror, R. O.; Dirks, R. M.; Grossman, J. P.; Xu, H.; Shaw, D. E. Biomolecular Simulation: A Computational Microscope for Molecular Biology. Annu. Rev. Biophys. 2012, 41, 429–452. (3) Lelièvre, T.; Rousset, M.; Stoltz, G. Free Energy Computations: A Mathematical Perspective; World Scientific: 2010. (4) Darve, E.; Pohorille, A. Calculating Free Energies Using Average Force. J. Chem. Phys. 2001, 115, 9169–9183. (5) Darve, E.; Wilson, M. A.; Pohorille, A. Calculating Free Energies Using a Scaled-Force Molecular Dynamics Algorithm. Mol. Simul. 2002, 28, 113–144. (6) Hénin, J.; Chipot, C. Overcoming Free Energy Barriers Using Unconstrained Molecular Dynamics Simulations. J. Chem. Phys. 2004, 121, 2904–2914. (7) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562–12566. (8) Grubmüller, H. Predicting Slow Structural Transitions in Macromolecular Systems: Conformational Flooding. Phys. Rev. E 1995, 52, 2893–2906. (9) Huber, T.; Torda, A. E.; Gunsteren, W. F. van. Local Elevation: A Method for Improving the Searching Properties of Molecular Dynamics Simulation. J. Comput.-Aided Mol. Des. 1994, 8, 695–708. (10) Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23, 187–199. (11) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011–1021. (12) Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. TRANSITION PATH SAMPLING: Throwing Ropes over Rough Mountain Passes, in the Dark. Annu. Rev. Phys. Chem. 2002, 53, 291–318. (13) Bolhuis, P. G.; Dellago, C.; Chandler, D. Reaction Coordinates of Biomolecular Isomerization. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 5877–5882. (14) Hazel, A.; Chipot, C.; Gumbart, J. C. Thermodynamics of Deca-Alanine Folding in Water. J. Chem. Theory Comput. 2014, 10, 2836–2844. (15) Geissler, P. L.; Dellago, C.; Chandler, D. Kinetic Pathways of Ion Pair Dissociation in Water. J.

29



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 30 of 36

Phys. Chem. B 1999, 103, 3706–3710. (16) Sugita, Y.; Okamoto, Y. Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 1999, 314, 141–151. (17) Sugita, Y.; Kitao, A.; Okamoto, Y. Multidimensional Replica-Exchange Method for Free-Energy Calculations. J. Chem. Phys. 2000, 113, 6042–6051. (18) Raiteri, P.; Laio, A.; Gervasio, F. L.; Micheletti, C.; Parrinello, M. Efficient Reconstruction of Complex Free Energy Landscapes by Multiple Walkers Metadynamics. J. Phys. Chem. B 2006, 110, 3533–3539. (19) Minoukadeh, K.; Chipot, C.; Lelièvre, T. Potential of Mean Force Calculations: A Multiple-Walker Adaptive Biasing Force Approach. J. Chem. Theory Comput. 2010, 6, 1008– 1017. (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, 5276–5285. (21) Wojtas-Niziurski, W.; Meng, Y.; Roux, B.; Bernèche, S. Self-Learning Adaptive Umbrella Sampling Method for the Determination of Free Energy Landscapes in Multiple Dimensions. J. Chem. Theory Comput. 2013, 9, 1885–1895. (22) Meng, Y.; Roux, B. Efficient Determination of Free Energy Landscapes in Multiple Dimensions from Biased Umbrella Sampling Simulations Using Linear Regression. J. Chem. Theory Comput. 2015, 11, 3523–3529. (23) Chen, M.; Yang, W. On-the-Path Random Walk Sampling for Efficient Optimization of Minimum Free-Energy Path. J. Comput. Chem. 2009, 30, 1649–1653. (24) 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. (25) Hamelberg, D.; Oliveira, C. A. F. de; McCammon, J. A. Sampling of Slow Diffusive Conformational Transitions with Accelerated Molecular Dynamics. J. Chem. Phys. 2007, 127, 155102. (26) Miao, Y.; Sinko, W.; Pierce, L.; Bucher, D.; Walker, R. C.; McCammon, J. A. Improved Reweighting of Accelerated Molecular Dynamics Simulations for Free Energy Calculation. J. Chem. Theory Comput. 2014, 10, 2677–2689.

30


Page 31 of 36

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


(27) Chipot, C.; Lelièvre, T. Enhanced Sampling of Multidimensional Free-Energy Landscapes Using Adaptive Biasing Forces. SIAM J. Appl. Math. 2011, 71, 1673–1695. (28) 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, 3506– 3513. (29) Lelièvre, T.; Rousset, M.; Stoltz, G. Computation of Free Energy Profiles with Parallel Adaptive Dynamics. J. Chem. Phys. 2007, 126, 134111. (30) Zheng, L.; Chen, M.; Yang, W. Random Walk in Orthogonal Space to Achieve Efficient Free-Energy Simulation of Complex Systems. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 20227– 20232. (31) Zheng, L.; Yang, W. Practically Efficient and Robust Free Energy Calculations: Double-Integration Orthogonal Space Tempering. J. Chem. Theory Comput. 2012, 8, 810–823. (32) Fiorin, G.; Klein, M. L.; Hénin, J. Using Collective Variables to Drive Molecular Dynamics Simulations. Mol. Phys. 2013, 111, 3345–3362. (33) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781–1802. (34) Comer, J.; Roux, B.; Chipot, C. Achieving Ergodic Sampling Using Replica-Exchange Free-Energy Calculations. Mol. Simul. 2014, 40, 218–228. (35) 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 2015, 119, 1129–1151. (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, 281–290. (37) Takekiyo, T.; Imai, T.; Kato, M.; Taniguchi, Y. Temperature and Pressure Effects on Conformational Equilibria of Alanine Dipeptide in Aqueous Solution. Biopolymers 2004, 73, 283– 290. (38) Grdadolnik, J.; Mohacek-Grosev, V.; Baldwin, R. L.; Avbelj, F. Populations of the Three Major Backbone Conformations in 19 Amino Acid Dipeptides. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 1794–1798.

31



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 32 of 36

(39) Balon, K.; Riebesehl, B. U.; Müller, B. W. Drug Liposome Partitioning as a Tool for the Prediction of Human Passive Intestinal Absorption. Pharm. Res. 1999, 16, 882–888. (40) Jämbeck, J. P. M.; Lyubartsev, A. P. Exploring the Free Energy Landscape of Solutes Embedded in Lipid Bilayers. J. Phys. Chem. Lett. 2013, 4, 1781–1787. (41) Plavec, J.; Tong, W.; Chattopadhyaya, J. How Do the Gauche and Anomeric Effects Drive the Pseudorotational Equilibrium of the Pentofuranose Moiety of Nucleosides? J. Am. Chem. Soc. 1993, 115, 9734–9746. (42) Huang, M.; Giese, T. J.; Lee, T.-S.; York, D. M. Improvement of DNA and RNA Sugar Pucker Profiles from Semiempirical Quantum Methods. J. Chem. Theory Comput. 2014, 10, 1538–1545. (43) Ensing, B.; Laio, A.; Parrinello, M.; Klein, M. L. A Recipe for the Computation of the Free Energy Barrier and the Lowest Free Energy Path of Concerted Reactions. J. Phys. Chem. B 2005, 109, 6676–6687. (44) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586–3616. (45) MacKerell, A. D.; Feig, M.; Brooks, C. L. Improved Treatment of the Protein Backbone in Empirical Force Fields. J. Am. Chem. Soc. 2004, 126, 698–699. (46) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM General Force Field: A Force Field for Drug-like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31, 671–690. (47) Yu, W.; He, X.; Vanommeslaeghe, K.; MacKerell, A. D. Extension of the CHARMM General Force Field to Sulfonyl-Containing Compounds and Its Utility in Biomolecular Simulations. J. Comput. Chem. 2012, 33, 2451–2468. (48) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33–38. (49) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method. J. Chem. Phys. 1995, 103, 4613–4621. (50) 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.

32


Page 33 of 36

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


(51) Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible Multiple Time Scale Molecular Dynamics. J. Chem. Phys. 1992, 97, 1990–2001. (52) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327–341. (53) Andersen, H. C. Rattle: A “velocity” Version of the Shake Algorithm for Molecular Dynamics Calculations. J. Comput. Phys. 1983, 52, 24–34. (54) Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952–962. (55) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N⋅log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089–10092. (56) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-4 Theory Using Reduced Order Perturbation Theory. J. Chem. Phys. 2007, 127, 124105. (57) Frisch, M.; Trucks, G.; Schlegel, H.; Scuseria, G.; Robb, M.; Cheeseman, J.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G.; et al. Gaussian 09, Gaussian. Inc., Wallingford, CT 2009. (58) Ling, S.; Gutowski, M. Different Conformations of 2′-Deoxycytidine in the Gas and Solid Phases: Competition between Intra- and Intermolecular Hydrogen Bonds. J. Phys. Chem. A 2016, 120, 8199-8210. (59) Lelièvre, T.; Minoukadeh, K. Long-Time Convergence of an Adaptive Biasing Force Method: The Bi-Channel Case. Arch. Rational Mech. Anal. 2011, 202, 1–34. (60) Lelièvre, T.; Rousset, M.; Stoltz, G. Long-Time Convergence of an Adaptive Biasing Force Method. Nonlinearity 2008, 21, 1155–1181. (61) Legoll, F.; Lelièvre, T. Effective Dynamics Using Conditional Expectations. Nonlinearity 2010, 23, 2131–2165. (62) Legoll, F.; Lelievre, T.; Olla, S. Pathwise Estimates for an Effective Dynamics. 2016. arXiv:1605.02644[math.PR]. arXiv.org e-Print archive. https://arxiv.org/abs/1605.02644 (accessed Dec. 7, 2016). (63) Alrachid, H.; Lelièvre, T. Long-Time Convergence of an Adaptive Biasing Force Method: Variance Reduction by Helmholtz Projection. SMAI J. Comput. Math. 2015, 1, 55–82.

33



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 34 of 36

(64) 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 [Online early

access].

DOI:

10.1021/acs.jpcb.6b10055.

Published

Online:

http://pubs.acs.org/doi/abs/10.1021/acs.jpcb.6b10055 (accessed Feb. 15, 2017).

34


Dec

13,

2016.

Page 35 of 36

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


TOC graphic:

35



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 36 of 36

The Extended Generalized Adaptive Biasing Force Algorithm for

Recommend Documents