Baixe Exploring Peptide Conformational Energy Landscapes: ABF & Metadynamics e outras Notas de estudo em PDF para Engenharia Elétrica, somente na Docsity!
Exploring Multidimensional Free Energy Landscapes
Using Time-Dependent Biases on Collective Variables
Je´rome He´nin,* ,†,|^ Giacomo Fiorin,‡^ Christophe Chipot,§,⊥^ and Michael L. Klein‡
Center for Molecular Modeling, Department of Chemistry, Uni V ersity of Pennsyl V ania, Philadelphia, Pennsyl V ania 19104, Institute for Computational Molecular Science and Department of Chemistry, Temple Uni V ersity, Philadelphia, Pennsyl V ania 19122, and Department of Physics and Beckman Institute for Ad V anced Science and Engineering, Uni V ersity of Illinois at Urbana - Champaign, Urbana, Illinois 61820 Received August 18, 2009
Abstract: A new implementation of the adaptive biasing force (ABF) method is described. This implementation supports a wide range of collective variables and can be applied to the computation of multidimensional energy profiles. It is provided to the community as part of a code that implements several analogous methods, including metadynamics. ABF and metady- namics have not previously been tested side by side on identical systems. Here, numerical tests are carried out on processes including conformational changes in model peptides and translocation of a halide ion across a lipid membrane through a peptide nanotube. On the basis of these examples, we discuss similarities and differences between the ABF and metadynamics schemes. Both approaches provide enhanced sampling and free energy profiles in quantitative agreement with each other in different applications. The method of choice depends on the dimension of the reaction coordinate space, the height of the barriers, and the relaxation times of degrees of freedom in the orthogonal space, which are not explicitly described by the chosen collective variables.
Introduction
A variety of approaches for accelerated sampling and mapping of free energy landscapes from molecular simula- tions have been proposed over the years (see refs 1–3 for reviews). Typically, these approaches have only been used by a limited number of groups that specialize in theory and method development. Only rarely have such methods been made readily available to the broad computational chemistry and biophysics community, as this requires well-documented implementations compatible with the standard tools of this
community: parallel simulation programs capable of high- throughput, large-scale calculations. An implementation of the metadynamics approach of Laio and Parrinello^4 for mainstream simulation packages has been made available only very recently.^5 The previous publicly available implementation^6 of the adaptive biasing force method7,8^ (ABF), in version 2.6 of NAMD,^9 has been applied successfully to a number of challenging cases. The domains of application of ABF include the recognition and association of peptides or proteins,10–12^ peptide- or protein-lipid interactions,13– small molecules interacting in a confined environment,^16 cyclodextrin association with cholesterol,^17 steroid drugs,^18 and molecular ions,^19 as well as cyclodextrin self-assembly.^20 Translocation of molecules or ions through natural, trans- membrane channel proteins21–24^ and transporters,^25 through synthetic pores,^26 and across simple liquid interfaces 27 have also been studied. Another class of applications involves conformational changes in peptides,28,29^ proteins,30–32^ and
- Corresponding author e-mail: [email protected]. † (^) University of Pennsylvania. ‡ (^) Temple University. § (^) University of Illinois at Urbana-Champaign. | (^) Current address: Laboratoire d’Inge´nierie des Syste`mes Mac-
romole´culaires, CNRS UPR9027, Marseille, France. ⊥ (^) On leave from Equipe de Dynamique des Assemblages
Membranaires, UMR CNRS/UHP 7565, Nancy Universite´ BP 239, Nancy, France.
J. Chem. Theory Comput. XXXX, xxx, 000 A
10.1021/ct9004432 XXXX American Chemical Society
nucleic acids.^33 Despite the number and variety of such applications, that implementation carried significant technical limitations, particularly its restriction to one-dimensional free energy profiles. When a low-dimension reduced representation comprised of a few degrees of freedom is used to describe a complex process, hidden barriers orthogonal to the chosen parameters are likely to exist. The orthogonal space random walk strategy34,35^ has been proposed by Zheng et al. as a means of overcoming such hidden barriers. While a promising idea, likely to be further explored and built upon in the future, it treats the orthogonal space using a single degree of freedom, which may or may not suffice to overcome hidden barriers effectively in complex examples; in many instances, a well- chosen degree of freedom may well yield better results. Still, the orthogonal space random walk approach has the advan- tage of generality, as it extends a predefined reaction coordinate space without requiring any additional physical insight into the particular process being examined.
Using collective variables well-adapted to the chemical or biophysical process under scrutiny is critical, and specific problems frequently require tailored variables. We recently developed a new code, the Collective Variables Module, for version 2.7 of the high-performance simulation program NAMD. This code supports a large set of commonly employed variables, offers the possibility to use polynomial combinations of such variables, and can be readily adapted to deal with atypical problems; the full list of features and technical details will be discussed elsewhere. In this contri- bution, we illustrate its most important application, that is, sampling multidimensional collective variable space and reconstructing free energy landscapes. Example simulations are presented, which make use of the first publicly released ABF implementation capable of multidimensional calcula- tions, and are discussed alongside the results of identical simulations performed with the metadynamics method.^4 In the following section, the theoretical framework un- derlying this ABF implementation is described, and its range of applicability as well as its technical limitations are discussed. Next, physical processes in four molecular systems are explored using ABF, conformational equilibria of N -acetyl- N ′-methylalaninamide (NANMA), Met-enkephalin, and deca- alanine, as well as ion diffusion through a membrane- spanning peptide nanotube. The metadynamics approach is also applied to the NANMA and nanotube examples. The deca-alanine case is used to document the application of ABF to a three-dimensional reaction coordinate. The choice of reaction coordinate space, numerical behavior, and conver- gence of the simulations, as well as compared properties of the two methods, are discussed.
Methods
Defining Reaction Coordinates. The strategy described here consists of using ABF or metadynamics to map a complex, slow molecular process, based on simulated trajectories that are orders of magnitude shorter than its natural time scale. This can be achieved by navigating a carefully chosen reduced representation, the “reaction co-
ordinate space”, in an accelerated fashion. The minimum requirement for this approach to be useful is that the reduced representation resolves the end points of the transformation and, more generally, all states that one wishes to describe based on empirical knowledge of the system. For numerical efficiency of the sampling scheme, however, the chosen degrees of freedom should capture all kinetically significant regions of configuration space: the metastable intermediates and most probable transition pathways. In chemical terminology, a reaction coordinate is a one- dimensional geometric parameter that can be used to measure the progression of a reaction.^36 Moving from the realm of chemical reactions to that of physical transformations in soft matter and biological systems, however, fluctuations along many degrees of freedom may become as important to the reaction kinetics as the progression along any particular pathway. One-dimensional descriptors then become less useful, while constructing single variables that may play this role becomes more cumbersome and less intuitive. For the purpose of numerical simulations, the optimal situation is to achieve time scale separation, whereby all key slow degrees of freedom are described explicitly, so that other degrees of freedom coupled to the transformation relax on a short time scale, as compared to the length of the simulated trajectories. This time scale influences both the diffusion rate of the system in reaction coordinate space and the rate of convergence of quantities that are measured as a function of the reaction coordinates, such as the free energy gradient in ABF calculations. Thermodynamic Integration in Configuration Space. This section gives a brief historical overview of the theoretical results that led to the ABF method. 7 The general principles of thermodynamic integration (TI) can be found in early work by Kirkwood 37 and Zwanzig. 38 In TI, the free energy derivative is computed as the ensemble average of an instantaneous force, F , acting on the reaction parameter : d A d
) -〈 F 〉 (1)
In most applications of TI to configurational variables, sampling along the reaction pathway is obtained by con- straining the reaction coordinate, the so-called “blue moon ensemble”. In one of the earliest cases, F was simply obtained as the force exerted by the solvent on two atomic ions, projected onto the interionic distance. 39 Shortly after- ward, a general expression for the average force was proposed by Carter et al.^40 The expression features a Jacobian correction term, purely geometric in origin, and is based on an explicit set of generalized coordinates (, q ) including the reaction coordinate :
F (, q ) ) - ∂ U (, q ) ∂
The explicit coordinate transform from ( xi ) to (, q ) is needed to define and compute both the Jacobian determinant | J (, q )| and the partial derivative ∂ U (, q )/∂ of the potential energy U. The latter quantity depends implicitly on the vector field (∂ xi /∂), hereafter referred to as “inverse gradient”. This vector field can be thought of as the direction along which
B J. Chem. Theory Comput., Vol. xxx, No. xx, XXXX He´nin et al.
could be recast into a more scalable form, paving the way for higher-dimension applications.
Computational Details
Molecular Dynamics Simulations. All simulations re- ported here were carried out using version 2.7b1 of the molecular dynamics program NAMD.9,47^ Condensed-phase simulations were performed in the isobaric-isothermal ensemble. The pressure and the temperature were fixed at 1 bar and 300 K, respectively, employing the Langevin piston algorithm^48 and softly damped Langevin dynamics. Periodic boundary conditions were applied in the three directions of Cartesian space. Short-range Lennard-Jones and Coulomb interactions were truncated smoothly by means of a 12 Å spherical cutoff with a switching function applied beyond 10 Å. The particle-mesh Ewald method^49 was employed to compute long-range electrostatic interactions. The Verlet I r -RESPA multiple time-step integrator^50 was used with a time step of 2 and 4 fs for for updating short- and long- range forces, respectively. Covalent bonds involving a hydrogen atom were constrained to their equilibrium length. Gas-phase simulations were performed using a 0.5 fs time step, which is appropriate to ensure energy conservation, and bond lengths were not constrained. Other parameters were similar to those of condensed-phase simulations. The dif- ferent chemical systems described in the present contribution were described by the all-atom CHARMM force field,^51 supplemented by the TIP3P water model. 52 Free Energy Calculations. The present results were obtained using a software framework known as Collective Variables Module and implemented in NAMD, versions 2.7b1 and following. Detailed user-oriented documentation is available;^47 technical details will be published elsewhere. Conformational Equilibrium of N-acetyl-N ′ -methylalanin- amide. The first application consists of a proof-of-concept simulation of the prototypical, terminally blocked amino acid N -acetyl- N ′-methylalaninamide (NANMA), often referred to as “alanine dipeptide”.^53 Conformational sampling was performed in vacuum. The φ and ψ torsional angles of the backbone were handled as coupled variables covering each the full, [- 180 °; + 180 °] range of the Ramachandran free energy map.^54 To increase the efficiency of the calculation, the latter map was split into four individual quadrants, corresponding to fully independent simulations. Each quad- rant was discretized into bins 2.5° × 2.5° wide, in which the force acting along the collective variables was accrued. In each quadrant, 25 ns of sampling was collected. A threshold of 100 force samples was set prior to application of the adaptive biasing force. Reconstruction of the complete free energy landscape was achieved by numerical integration of the two-dimensional gradients. The Ramachandran map was also sampled by means of the metadynamics algorithm. Sampling was collected from a 30 ns trajectory, and Gaussian biasing potentials of width 10° and height 0.1 kcal/mol were accumulated every 500 fs. The free energy difference between the C 7 eq and C 7 ax states was computed by integration over the corresponding regions V and V ′ of (φ, ψ)-space according to:
e - ∆ A^ )
∫ V e
∫ V' e
Transition between the two lowest free energy states of the Ramachandran map, that is, C (^7) eq and C (^7) ax , which are stabilized by an intramolecular hydrogen bond formed between the carbonyl moiety of one terminus and the amino moiety of the other, was investigated with one-dimensional ABF, using as a collective variable the difference between two distance root mean-square deviations, ) rmsd(C (^7) eq )
- rmsd(C 7 ax ). To ensure orthogonality of the variables according to eq 6, the two RMSDs were defined as two distinct subsets formed of six atoms of the peptide chain. Three independent, 20, 20, and 40 ns long, simulations were run, using a threshold of 5000 force samples prior to applying the adaptive biasing force along the chosen degrees of freedom. Chloride Ion Permeation across a Peptide Nanotube. In this second application, translocation of an halide ion through a chemically engineered tubular structure is examined. Peptide nanotubes, which can be viewed as tailored synthetic ion channels, result from the self-assembly of cyclic peptides formed by alternated D - L - R-amino acids,55,56^ by means of a network of intermolecular hydrogen bonds.^57 The peptide nanotube considered here consisted of eight stacked cyclo [LW] 4 units, where underlined letters denote D-amino acids, immersed in a thermalized palmitoyl-oleoyl- phosphatidylcholine (POPC) bilayer formed by 48 lipid units, in equilibrium with 1572 water molecules. The complete molecular assembly was placed in a simulation cell of initial dimensions equal to 36 × 41 × 79 Å^3. The two-dimensional free energy landscape delineating the translocation of a chloride ion across the tubular structure was determined along the longitudinal, , and the radial, F, directions of the latter. Specifically, the model reaction coordinate was chosen as a subset of cylindrical polar coordinates: the distance separating the ion from the center of mass of the peptide nanotube, projected onto its long axis, in conjunction with the distance separating the ion from the axis. The reaction path spanned 40 and 3 Å in the - and in the F-directions, respectively. In the ABF calculation, force samples were accrued in bins 0. Å wide. To increase the efficiency of the calculation, it was stratified into four nonoverlapping windows extending over 10 Å each in the -direction and in which individual 30 ns trajectories were generated, corresponding to a total simula- tion time of 120 ns. A metadynamics simulation was performed using the same pair of variables, similarly split into four windows along z. Gaussian hills of width 0.3 Å and height 0.1 kcal/mol were added every 200 fs; the calculation was run for 44 ns. Structure of Met-enkephalin in an Aqueous Solution. In this third application, a set of collective variables is utilized to explore the possible conformations of the short peptide Met-enkephalin in an explicit solvent. Met-enkephalin is an endogenous opioid, five-residue neurotransmitter peptide, YGGFM, found in mammals and known to inhibit the release of neurotransmitters upon activation of the relevant opioid receptors. On account of its small size and biological relevance, it has served as a paradigmatic system for
D J. Chem. Theory Comput., Vol. xxx, No. xx, XXXX He´nin et al.
conformational search based on a variety of computational approaches.58–
The molecular system consisted of Met-enkephalin im- mersed in a bath of 778 water molecules, which corresponds to a simulation cell of initial dimensions equal to 29 × 29 × 29 Å^3. Conformational search was conducted in a two- dimensional space, combining the radius of gyration ( R g) of the short peptide to its distance rmsd with respect to a reference, helical structure. The reaction path spanned, respectively, 3.5 and 4.0 Å, in the R g and rmsd directions. To ensure that the force measured along one variable does not act on the other (eq 6), two distinct subsets of atoms were defined to compute the distance rmsd and the radius of gyration, the five R-carbon atoms and all other heavy atoms of the peptide chain, respectively. The instantaneous force was accrued in bins 0.05 Å wide. No adaptive biasing force was applied below a threshold of 200 samples. The limited range covered by the two variables obviated the need for a stratification strategy. The two-dimensional free energy landscape reported in the present contribution was obtained from a total simulation time of 80 ns. Conformational Landscape of Deca-alanine. To explore metastable conformations of deca-alanine in vacuum, starting points were chosen manually from preliminary ABF trajec- tories, and used as the seed for unbiased, 30 ns simulations, where relaxation and possible transitions to other local minima were monitored. Conformations typical of the long- lived, metastable conformers were extracted from these unbiased trajectories. The conformational free energy land- scape of the peptide was then explored by means of ABF- biased simulations in one, two, and three dimensions. The one-dimensional calculation was based on the end-to-end distance d , that is, the distance between the carbonyl carbon atoms of residues 1 and 10. The duration of the simulated trajectory was 500 ns. The two-dimensional, 200 ns ABF calculation was based on both d and the minimal rmsd between the current structure and a typical -hairpin conformation (structure E in Figure 6). The rmsd was calculated using R-carbon atoms only. Finally, a 400 ns three- dimensional ABF simulation was carried out: the set of collective variables was composed of three RMSDs, with respect to the typical R-helical, turn/3 10 -helix, and ω-shaped turn, respectively. The structures of these conformers are described in the Results section. The three RMSDs were defined on the basis of R-carbon, carbonyl carbon, and peptide nitrogen atoms, respectively; these nonoverlapping sets of atoms ensured that eq 6 was obeyed.
Results and Discussion
Conformational Equilibrium of N -Acetyl- N ′ -methyl- alaninamide. In the past 30 years, the conformational equilibrium of NANMA has been investigated in ample detail, employing a variety of numerical approaches and potential energy functions.4,8,53,73–93^ Accurate reproduction of the gas-phase, two-dimensional Ramachandran free energy map, therefore, hardly constitutes a methodological prowess. Popular, well-established approaches like umbrella sam- pling^94 used in conjunction with the weighted histogram
analysis method^95 have proven to be perfectly adapted to rank with the desired accuracy the conformational states accessible to this paradigmatic peptide.^82 Furthermore, in several instances, NANMA has served as a discriminating test case for assessing the performance of novel numerical schemes.8,83, In the present work, determination of the two-dimensional free energy landscape of NANMA only represents the necessary preamble to an independent series of computations based on variables of higher collectivity than the mere φ and ψ torsional angles of the backbone. However predictable, the results of these preliminary free energy calculations depicted in Figure 1 agree quantitatively with previous investigations; see, for instance, ref 96. In particular, (φ, ψ)- integration over the basins corresponding to the C 7 eq and C 7 ax conformations, characterized by a long-lived intramolecular hydrogen bond, yields a free energy difference of 2.5 kcal/ mol, in favor of the equatorial motif.78, The metadynamics calculation yields a free energy land- scape identical to that determined by means of ABF (data not shown). Fast-relaxing, low-friction systems such as NANMA do not place high demands on the algorithms used to sample them, as demonstrated by early successes obtained
Figure 1. (A) Two-dimensional (φ, ψ) free energy landscape of N -acetyl- N ′-methylalanylamide in the gas phase derived from a reference ABF calculation. (B) Free energy change along collective variable ) rmsd(C (^7) eq ) - rmsd(C (^7) ax ), the difference between the distance rmsd with respect to the C 7 eq and C 7 ax conformational states. The free energy profiles are obtained from independent 20 (light curve), 20 (dark curve), and 40 ns (black curve) ABF simulations.
Multidimensional Free Energy Landscapes J. Chem. Theory Comput., Vol. xxx, No. xx, XXXX E
midpoint between two adjacent cyclic peptides, the separation of which amounts to about 4.7 Å. Congruent with the simulations of the sodium ion,^26 in-plane coordination appears to be also the energetically preferred state for chloride. It is remarkable that the halide ion does not hop between cyclic peptides following a rectilinear path, collinear to the long axis of the synthetic channel at ) 0 Å, but instead grazes the wall of the latter to form intermittent interactions with peptide amino moieties. This diffusion mode is reflected in a curvilinear low-free energy path apparent in Figure 2. Structure of Met-enkephalin in an Aqueous Solution. On account of its unusual flexibility, Met-enkephalin can bind interchangeably either the δ, the κ, or the μ-opioid receptor, where it is anticipated to adopt distinct conformations. This, in large measure, rationalizes why experimental approaches have hitherto failed to propose a converging view of the native structure of this short peptide.26,106,107^ Although it is reasonable to believe that the marked flexibility of the backbone would be reflected in a generally flat free energy landscape, equilibrium MD simulations at a single temper- ature have proven to prevent exhaustive exploration of conformational space, the peptide chain being recurrently trapped kinetically in a random, unrepresentative conforma- tion.^60 Despite its reduced length, Met-enkephalin, therefore, constitutes a pathological case for standard MD, but a well- suited candidate for more recent, multicanonical approaches, like replica-exchange MD,^108 aimed at enhancing ergodic sampling. The two-dimensional free energy map of Figure 3 sheds light on the conformational space accessible to the short peptide along the directions of the position rmsd with respect to an arbitrary chosen helical motif and of the radius of gyration of the backbone. As has been commented on previously,^67 the free energy landscape of Met-enkephalin in aqueous solution consists of a rather wide, shallow basin, featuring a number of pronounced minima. At low distance rmsd and R g, a compact, 3 10 -helix motif ( I ) is observed, resulting primarily from the formation of an intramolecular hydrogen bond between the carbonyl group of residue i and
the amino group of residue i + 3, albeit transitory intercon- version to an R-helical i to i + 4 hydrogen bond can be detected. Still at low R g, but at a larger distance rmsd, around 2.8 Å, a free energy minimum ( II ) corresponding to an embryonic helix turn emerges about 2.1 kcal/mol higher than that of conformation I. A third minimum ( III ) can be found at a distance rmsd of 2.0 Å and an R g of 5.2 Å, with a free energy equal to that of the helical motif I within statistical accuracy. Its U-shaped structure is essentially nonhelical, featuring an intermittent hydrogen bond formed between the terminal, blocking amino moiety and the carbonyl group of the first glycine residue. Last, a fourth conformational minimum ( IV ) can be found at a distance rmsd of 1.4 Å and an R g of 4.8 Å, and corresponds to a free energy only 0.6 kcal/mol above that of structures I and III. In this free energy minimum, the peptide chain adopts a γ-turn confor- mation, often encountered in -hairpin motifs. It is apparent from the present results that, using a reduced set of collective variables and a single temperature, ABF is able to recover the complete free energy landscape of Met- enkephalin, virtually identical to that reported by Sanbon- matsu and Garcia on the basis of parallel-tempering MD simulations.^67 The (rmsd, R g) two-dimensional map confirms, indeed, that this free energy landscape is relatively shallow and consists of essentially four distinct ensembles of conformations characterized by overall comparable free energies, barring structure II , and separated by appreciably low barriers. The local minima are primarily distinguished by their rmsd value. Therefore, in this case, the second collective variable, R g, is not essential to resolve the low free energy states, but it accelerates the sampling of transition pathways between them, and therefore the overall convergence of the calcula- tion. The rmsd could be considered a one-dimensional reaction coordinate, while R g is a degree of freedom in the orthogonal space that is included in the ABF calculation for numerical efficiency. Conformational Free Energy Landscape of Deca-ala- nine. Deca-alanine in vacuum has been used several times as a computationally inexpensive and seemingly simple toy
Figure 3. Two-dimensional (rmsd, R g ) free energy landscape of Met-enkephalin in aqueous solution. The position rmsd with respect to a helical conformation is measured over R-carbon atoms only. The radius of gyration, R g, is evaluated over all heavy, but R-carbon atoms. Contours of the main two-dimensional map are separated by 0.5 kcal/mol. The I - IV key conformations of Met-enkephalin are charted on the inset free energy landscape and depicted to the right as cartoon and ball-and-stick representations. In the cartoon representation, color indicates secondary structure. Molecular graphics were rendered with VMD.^109
Multidimensional Free Energy Landscapes J. Chem. Theory Comput., Vol. xxx, No. xx, XXXX G
model for conformational free energy calculations on peptides.6,28,110–113^ Stretching the peptide from its native, R-helical state to extended conformations, and back, is indeed a fairly straightforward transformation, for which relatively few pathways are possible. All simulations cited above, using the end-to-end distance of the peptide chain as a coordinate and various biasing schemes, achieve convergence and find potentials of mean force in good agreement with one another. Such simulations, however, only explore a small fraction of the conformational space available to deca-alanine. Indeed, when an attempt was made to sample the more compact conformations corresponding to end-to-end distances shorter than the R-helix length (i.e., less than 12 Å), it was found that several local minima, and the complex pathways linking them, were not resolved by that simple distance parameter. As a result, the complete free energy landscape of the peptide could not be mapped.^28 As can be seen from Figure 4B, using ABF to sample the end-to-end coordinate d provides very uniform histograms for elongated structures, and uneven, sporadic sampling of the compact region. Thus, the potential of mean force shown in Figure 4A, while very accurate for d greater than 12 Å, is
not converged for smaller values. Furthermore, it does not exhibit clear features that could be related to the metastable conformations that can be identified in that region. The reason behind such nonuniform sampling of the low- d range is that d is a highly degenerate descriptor of most of the conformational space of deca-alanine. Even with an ideally converged adaptive bias canceling the free energy profile, some of the transitions would remain rare events, essentially because the enthalpic or entropic barriers are orthogonal to the direction of the bias. Unbiased simulations starting from structures chosen from the ABF simulations indicate several metastable conformers whose lifetime is at least 30 ns. In order of increasing relative potential energy, these structures are the R-helix (0 kcal/ mol), a loop containing a short 3 10 -helical segment (3.5 kcal/ mol), an ω-shaped loop composed of and other turns (6. kcal/mol), a -bridged loop (9.3 kcal/mol), a -hairpin (9. kcal/mol), and a slightly expanded π-helix featuring transient i to i + 5 hydrogen bonds (11 kcal/mol). These conforma- tions are represented in Figure 6. To resolve the compact states that are merged in the one- dimensional landscape, the distance parameter was supple- mented with a second degree of freedom, the position rmsd from a typical -hairpin conformation. The results are represented in Figure 5. Although a large fraction of the reduced space is visited, sampling tends to accumulate in a few localized regions (Figure 5B), suggesting a failure of the collective variables to describe all of the free energy barriers in these regions. The global free energy minimum corresponding to the R-helix is clearly visible in the measured free energy landscape (Figure 5A). In contrast, the basins corresponding to metastable structures are almost coalesced and separated by low barriers, inconsistent with the long lifetimes of these conformers. Despite its limits, the two- dimensional profile is consistent with the one-dimensional free energy of 4 in the region containing the R-helix and the extended states. Finally, a set of three RMSDs was chosen to try and resolve the numerous conformers of similar compactness of deca-alanine. The three-dimensional PMF is shown in Figure
- In general, the differences in free energy between the metastable conformers and the R-helix are smaller than the potential energy differences, indicating significant entropy- enthalpy compensation. Indeed, the R-helix is more rigid than the other conformers. Because of the width of basins corresponding to individual conformers, and the short distance separating them in the reduced space, some con- formers are still not perfectly resolved, in particular in the low free energy region containing the -hairpin and -bridge structures. Indeed, the local free energy minimum in that region contains conformations similar to both structures, whose similarity renders difficult to distinguish in a low- dimension representation. The free energy landscape is, however, sufficiently resolved for the complete conforma- tional space to be explored at a much increased rate due to the adaptive bias. Here, because the full conformation space of flexible peptide is explored, the rmsd coordinate is not robust over its entire range. Specifically, at very high rmsd values, the structure becomes too dissimilar from the
Figure 4. (A) One-dimensional free energy profile of deca- alanine in vacuum, as a function of the end-to-end distance. Intermediate results are plotted for every 20 ns period as gray lines; the final result of the 500 ns ABF simulation is shown as a dark line. (B) Sampling histograms (on a logarithmic scale) corresponding to the one-dimensional ABF simulation of deca-alanine. Each line is the full histogram at the end of a 20 ns period of the simulation.
H J. Chem. Theory Comput., Vol. xxx, No. xx, XXXX He´nin et al.
observed effect of slowly relaxing hidden degrees of freedom is a change of the sampling behavior as the simulation progresses. Over time, the relative weight of newly accrued data in the average decreases, as does the responsiveness of the time-dependent bias to changes in the locally measured thermodynamic force. This is the desired behavior in the ideal context of a perfect reaction coordinate space, where all orthogonal degrees of freedom are either fast, or kinetically trapped and not relevant to the transformation under study. In real cases, the symptom is an initially efficient exploration of reduced space, followed by an increasing tendency to get trapped in local minima as they are “discovered” through relaxation in the orthogonal space and crossing of hidden barriers. This is observed, to a minor extent, in the 2D and 3D deca-alanine simulations described in the previous sections. Arguably, such behavior is indicative of failure of the reduced space to capture the reactive intermediates and pathways to a degree that allows for sampling over the intended (or technically feasible) time scale. Metadynamics is designed as a nonequilibrium method, whereas ABF relies on equilibrium sampling from the canonical ensemble. In practice, ABF simulations go through an initial stage during which the estimate of the local free energy gradient evolves rapidly; then the running average is updated on a much slower scale, and eventually stabilizes altogether. Applying a biasing force based on the early, fast- fluctuating estimate of the free energy gradient may push
the dynamics into a nonequilibrium regime. This initial departure from equilibrium conditions can be greatly reduced in practice by delaying the introduction of the biasing force until after the large fluctuations of the running average have ended: this is accomplished by waiting for a preset number of samples to be accrued in a local bin. 6,8,28^ Should any nonequilibrium effects occur, however, force samples mea- sured in that phase are still taken into account in the average, and, as mentioned above, this contribution never vanishes. Any bias (typically, an overestimated free energy gradient) due to irreversible work performed in the beginning of the simulation will taint the final data, with a weight that only decays as 1/ t. In contrast, exponential convergence of ABF has only been proved under somewhat restrictive assump- tions, including the use of a large number of replicas. 116 When a metadynamics run is carried out with constant hill height, irreversible work is performed throughout the simula- tion. This contribution may only be deemed to have been eliminated from the measured PMF once multiple sweeps have been observed (or, in higher dimension, multiple transitions between neighboring basins). Both issues (slow hidden variables and irreversible work) can be alleviated by restricting the data set used to reconstruct the free energy, eliminating selected, problematic data. In metadynamics, this is achieved by stopping the simulation once the system leaves the relevant region of configuration space. In ABF, data collected in the initial, out-of-equilibrium
Figure 6. Three-dimensional conformational free energy landscape of deca-alanine, as a function of the rmsd with respect to ideal R-helical, ω-shaped γ-turn, and turn/3 10 -helix conformations, respectively. The free energy was computed by numerically integrating the gradient obtained from a 400 ns ABF simulation. Free energy isosurfaces are shown for respective values of 3 (black), 6 (blue), 9 (green), 12 (yellow), and 24 kcal/mol (white). Average positions corresponding to metastable conformations are indicated by red spheres. Conformers are shown as ball-and-stick and cartoon representations, with color indicating secondary structure (purple, R-helix; blue, 3 10 -helix; red, π-helix; gray, turn; yellow, -strand; tan, -bridge; white, unstructured). Hydrogen bonds are indicated as blue lines. The conformers are labeled using red capital letters: (A) R-helix, (B) turn/3 10 -helix, (C) ω-shaped turn, (D) -bridge, (E) -hairpin, (F) π-helix, and (G) extended form (unstable).
J J. Chem. Theory Comput., Vol. xxx, No. xx, XXXX He´nin et al.
stage can be removed from the final average. Metadynamics performs well in low-friction systems such as the peptide nanotube example, because the filling rate (hill height divided by the hill accumulation period) can be increased while maintaining a level of irreversible work that does not hinder convergence.
Conclusion
An implementation of the ABF approach supporting multi- dimensional reaction coordinates based on sophisticated variables is proposed. This implementation is available in the scalable molecular-dynamics program NAMD, but can be readily incorporated into any MD platform. The method and its implementation have been tested on a variety of biomolecular systems. Multidimensional reaction descriptors improve the level of detail in which molecular processes can be mapped, as in the case of chloride permeation through a self-assembled peptide nanotube. In systems featuring par- ticularly complex free energy landscapes, such as the multiple metastable conformations of deca-alanine, efficient sampling is only possible when applying the adaptive bias in a two- or three-dimensional reduced space. The work of constructing an appropriate reduced representation is made easier by the availability of variables such as the position rmsd, and the flexibility afforded by linear combinations of predefined variables.
Two of the test cases were used to compare directly ABF and the metadynamics method, applied on the same variables, under otherwise identical simulation conditions. Both ap- proaches yield the expected results in terms of phase space exploration and sampling, and reconstruction of the free energy landscape. Neither can be said to be more efficient or accurate than the other on general grounds. Still, meta- dynamics does offer an inexpensive way to rapidly explore simple, fast-relaxing systems that are robust enough to withstand a high energy input, as is apparent from the toy model of NANMA, while ABF may constitute a safer option for more fragile systems. Differences are likely to become most apparent in difficult scenarios, where the two ap- proaches react differently to incomplete sampling and nonequilibrium effects. Whereas metadynamics is much more likely to perform significant irreversible work in late stages of such simulations, ABF may become inefficient and sample reduced regions for disproportionate times.
The most crucial issue remains the choice of a reaction coordinate space. As illustrated in the deca-alanine example, describing a small system with seemingly limited phase space may require multiple parameters, fine-tuned to resolve nearby points of conformation space that, nevertheless, belong to different low free energy basins or pathways. More sophis- ticated collective variables will be developed to describe complex, frustrated free energy landscapes. The existing algorithms will evolve toward better scalability, or be replaced altogether by new sampling methods. In any event, it is essential for the progress of large-scale applications that new developments in sampling and free energy algorithms be kept in phase with state-of-the-art, parallel simulation software.
Acknowledgment. J.H., G.F., and M.L.K. acknowledge funding from the National Institutes of Health. We are grateful to Dr. Tony Lelièvre for enlightening discussions, and to Drs. Christopher Harrison and Peter Freddolino for their help with the implementation and testing of numerical methods. References (1) Rodinger, T.; Pome`s, R. Curr. Opin. Struct. Biol. 2005 , 15 , 164–170. (2) Free Energy Calculations. Theory and Applications in Chemistry and Biology ; Chipot, C., Pohorille, A., Eds.; Springer Verlag: New York, 2006. (3) Lei, H.; Duan, Y. Curr. Opin. Struct. Biol. 2007 , 17 , 187–
(4) Laio, A.; Parrinello, M. Proc. Natl. Acad. Sci. U.S.A. 2002 , 99 , 12562–12565. (5) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: a portable plugin for free-energy calculations with molecular dy- namics; 2009; arXiv:0902.0874v3 [physics.comp-ph]; arX- iv.org ePrint archive, http://www.arXiv.org/abs/0902. (accessed Nov. 9, 2009). (6) He´nin, J.; Chipot, C. J. Chem. Phys. 2004 , 121 , 2904–2914. (7) Darve, E.; Pohorille, A. J. Chem. Phys. 2001 , 115 , 9169–
(8) Darve, E.; Rodrı´guez-Go´mez, D.; Pohorille, A. J. Chem. Phys. 2008 , 128 , 144120. (9) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhor- shid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale´, L.; Schulten, K. J. Comput. Chem. 2005 , 26 , 1781–1802. (10) He´nin, J.; Pohorille, A.; Chipot, C. J. Am. Chem. Soc. 2005 , 127 , 8478–8484. (11) Xu, J.; Crowley, M. F.; Smith, J. C. Protein Sci. 2009 , 18 , 949–959. (12) Halling-Brown, M. D.; Moss, D. S.; Sansom, C. E.; Shepherd, A. J. Philos. Trans. R. Soc. London, Ser. A 2009 , 367 , 2705–2716. (13) Weronski, P.; Jiang, Y.; Rasmussen, S. Biophys. J. 2007 , 92 , 3081–3091. (14) Gorfe, A. A.; Babakhani, A.; McCammon, J. A. Angew. Chem., Int. Ed. 2007 , 46 , 8234–8237. (15) Gorfe, A. A.; McCammon, J. A. J. Am. Chem. Soc. 2008 , 130 , 12624–12625. (16) Vaitheeswaran, S.; Thirumalai, D. Proc. Natl. Acad. Sci. U.S.A. 2008 , 105 , 17636–17641. (17) Yu, Y.; Chipot, C.; Cai, W.; Shao, X. J. Phys. Chem. B 2006 , 110 , 6372–6378. (18) Cai, W.; Sun, T.; Liu, P.; Chipot, C.; Shao, X. J. Phys. Chem. B 2009 , 113 , 7836–7843. (19) Rodriguez, J.; Elola, M. D. J. Phys. Chem. B 2009 , 113 , 1423–1428. (20) Rodriguez, J.; Semino, R.; Laria, D. J. Phys. Chem. B 2009 , 113 , 1241–1244. (21) Treptow, W.; Tarek, M. Biophys. J. 2006 , 90 , L64–L66. (22) Lamoureux, G.; Klein, M. L.; Berneche, S. Biophys. J. 2007 , 92 , L82–L84.
Multidimensional Free Energy Landscapes J. Chem. Theory Comput., Vol. xxx, No. xx, XXXX K
(79) Pellegrini, M.; Grønbech-Jensen, N.; Doniach, S. J. Chem. Phys. 1996 , 104 , 8639. (80) Neria, E.; Fischer, S.; Karplus, M. J. Chem. Phys. 1996 , 105 , 1902. (81) Bartels, C.; Karplus, M. J. Comput. Chem. 1997 , 18 , 1450–
(82) Chipot, C.; Pohorille, A. J. Phys. Chem. B 1998 , 102 , 281–
(83) Smith, P. E. J. Chem. Phys. 1999 , 111 , 5568–5579. (84) Apostolakis, J.; Ferrara, P.; Caflisch, A. J. Chem. Phys. 1999 , 110 , 2099–2108. (85) Bolhuis, P. G.; Dellago, C.; Chandler, D. Proc. Natl. Acad. Sci. U.S.A. 2000 , 97 , 5877–5882. (86) Andricioaei, I.; Dinner, A. R.; Karplus, M. J. Chem. Phys. 2003 , 118 , 1074–1084. (87) Chekmarev, D.; Ishida, T.; Levy, R. M. J. Phys. Chem. B 2004 , 108 , 19487–19495. (88) Jang, H.; Woolf, T. B. J. Comput. Chem. 2006 , 27 , 1136–
(89) Branduardi, D.; Gervasio, F. L.; Parrinello, M. J. Chem. Phys. 2007 , 126 , 054103. (90) Neale, C.; Rodinger, T.; Pomes, R. Chem. Phys. Lett. 2008 , 460 , 375–381. (91) Kwac, K.; Lee, K. K.; Han, J. B.; Oh, K. I.; Cho, M. J. Chem. Phys. 2008 , 128 , 105106. (92) Feig, M. J. Chem. Theory Comput. 2008 , 4 , 1555–1564. (93) Velez-Vega, C.; Borrero, E. E.; Escobedo, F. A. J. Chem. Phys. 2009 , 130 , 225101. (94) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977 , 23 , 187–199. (95) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. J. Comput. Chem. 1992 , 13 , 1011–1021. (96) Rosso, L.; Abrams, J. B.; Tuckerman, M. E. J. Phys. Chem. B 2005 , 109 , 4162–4167. (97) Sa´nchez-Quesada, J.; Isler, M. P.; Ghadiri, M. R. J. Am. Chem. Soc. 2002 , 124 , 10004–10005. (98) Kim, H. S.; Hartgerink, J. D.; Ghadiri, M. R. J. Am. Chem. Soc. 1998 , 120 , 4417–4424.
(99) Tarek, M.; Maigret, B.; Chipot, C. Biophys. J. 2003 , 85 , 2287–2298. (100) Asthagiri, D.; Bashford, D. Biophys. J. 2002 , 82 , 1176–
(101) Hwang, H.; Schatz, G. C.; Ratner, M. A. J. Phys. Chem. B 2006 , 110 , 6999–7008. (102) Hwang, H.; Schatz, G. C.; Ratner, M. A. J. Phys. Chem. B 2006 , 110 , 26448–26460. (103) Mamonov, A. B.; Coalson, R. D.; Nitzan, A.; Kurnikova, M. G. Biophys. J. 2003 , 84 , 3646–3661. (104) Aguilella-Arzo, M.; Aguilella, V. M.; Eisenberg, R. S. Eur. Biophys. J. 2005 , 34 , 314–322. (105) Marcus, Y. Ion Sol V ation ; Wiley and Sons: London, 1985.
(106) Khaled, M. A.; Long, M. M.; Thompson, W. D.; Bradley, R. J.; Brown, G. B.; Urry, D. W. Biochem. Biophys. Res. Commun. 1976 , 76 , 224–231. (107) Smith, G. D.; Griffin, J. F. Science 1978 , 199 , 1214–1216. (108) Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 1999 , 314 , 141–
(109) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996 , 14 , 33–8, 27-8. (110) Kumar, S.; Payne, P. W.; Va´squez, M. J. Comput. Chem. 1995 , 17 , 1269–1275.
(111) Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. J. Chem. Phys. 2003 , 119 , 3559–3566. (112) Minh, D. D. L.; McCammon, J. A. J. Phys. Chem. B 2008 , 112 , 5892–5897. (113) Forney, M. W.; Janosi, L.; Kosztin, I. Phys. Re V_. E_ 2008 , 78 , 051913. (114) Laio, A.; Rodriguez-Fortea, A.; Gervasio, F. L.; Ceccarelli, M.; Parrinello, M. J. Phys. Chem. B 2005 , 109 , 6714–6721. (115) Bussi, G.; Laio, A.; Parrinello, M. Phys. Re V_. Lett._ 2006 , 96 , 090601. (116) Lelièvre, T.; Rousset, M.; Stoltz, G. J. Chem. Phys. 2007 , 126 , 134111. CT
Multidimensional Free Energy Landscapes J. Chem. Theory Comput., Vol. xxx, No. xx, XXXX M
