ACS Publications. Most Trusted. Most Cited. Most Read
My Activity
CONTENT TYPES

Figure 1Loading Img

Data Reweighting in Metadynamics Simulations

  • Timo M. Schäfer
    Timo M. Schäfer
    Institut für Physik, Johannes Gutenberg-Universität Mainz, Mainz, Germany
    Graduate School Materials Science in Mainz, Mainz, Germany
  •  and 
  • Giovanni Settanni*
    Giovanni Settanni
    Institut für Physik, Johannes Gutenberg-Universität Mainz, Mainz, Germany
    Max Planck Graduate Center mit der Johannes Gutenberg-Universität Mainz, Mainz, Germany
    *Email: [email protected]
Cite this: J. Chem. Theory Comput. 2020, 16, 4, 2042–2052
Publication Date (Web):March 19, 2020
https://doi.org/10.1021/acs.jctc.9b00867

Copyright © 2020 American Chemical Society. This publication is licensed under these Terms of Use.

  • Open Access

Article Views

4946

Altmetric

-

Citations

LEARN ABOUT THESE METRICS
PDF (5 MB)
Supporting Info (1)»

Abstract

The data collected along a metadynamics simulation can be used to recover information about the underlying unbiased system by means of a reweighting procedure. Here, we analyze the behavior of several reweighting techniques in terms of the quality of the reconstruction of the underlying unbiased free energy landscape in the early stages of the simulation and propose a simple reweighting scheme that we relate to the other techniques. We then show that the free energy landscape reconstructed from reweighted data can be more accurate than the negative bias potential depending on the reweighting technique, the stage of the simulation, and the adoption of well-tempered or standard metadynamics. While none of the tested reweighting techniques from the literature provides the most accurate results in all the analyzed situations, the one proposed here, in addition to helping simplifying the reweighting procedure, converges quickly and precisely to the underlying free energy surface in all the considered cases, thus allowing for an efficient use of limited simulation data.

1. Introduction

ARTICLE SECTIONS
Jump To

Molecular dynamics simulations allow following the behavior of matter at the molecular level with subnanometer resolution, thus representing an important tool for the dissection of physical phenomena in many different fields. Their applicability is limited by the accuracy of the underlying force fields, which are only an approximation of the real physical systems, and by the size and time scales of the phenomena under scrutiny, which may be too large for the available computer resources. (1,2) To address the latter point, besides advances in the available hardware and software, (2) many different techniques (3) have been developed to improve the sampling efficiency including, but in no way limited to, umbrella sampling, (4) replica exchange molecular dynamics, (5−7) path sampling techniques, (8) adaptive biasing force, (9) temperature-accelerated molecular dynamics, (10) and metadynamics. (11−13) In particular, metadynamics allows for an estimate of the free energy landscape of large multidimensional systems as a function of a small number of relevant collective variables (CV). According to this technique, during the course of the simulation, a time-dependent bias potential is built up in the space of the CV as a sum of repulsive Gaussians centered on the points visited along the trajectory. The accumulation of the bias potential facilitates a thorough exploration of the phase space allowing for the crossing of high free energy barriers. After convergence, the free energy of the system as a function of the CV can be recovered from the bias potential built during the simulation. (11,13) The conformations of the system sampled during the metadynamics simulation can be used to obtain information about the unbiased state of the system, including unbiased estimate of observables, by means of a reweighting procedure. Data reweighting can also be used to assess the free energy of the system as a function of the biased variables in alternative to the bias potential, (14−16) and this represents a common way to test the accuracy of the reweighting.
Several reweighting procedures have been proposed in the literature. (14−18) In particular, the reweighting schemes proposed by Tiwary and Parrinello (16) and by Branduardi et al. (15) provide a simple and general implementation included in the program PLUMED, (19,20) a popular free-energy and metadynamics simulation package, and described in its tutorial courses. The method by Bonomi et al. (14) has also been made available in an earlier version of PLUMED. While the methods by Branduardi et al. (15) and Bonomi et al. (14) are specifically indicated for well-tempered metadynamics, which is a variant of metadynamics where the height of the repulsive Gaussians is gradually reduced along the simulation, (13) the one from Tiwary and Parrinello (16) is also indicated for standard metadynamics. The assumptions behind the three reweighting schemes make them effective after a transient period from the beginning of the simulation. Recently, several tests have been conducted on these reweighting schemes focusing on their behavior on the long time scales. (21) Often, however, in real case scenarios involving large and complex molecular systems, we are faced with the problem of having a limited amount of trajectory data. Thus, it is important to evaluate how quickly each reweighting scheme allows for reaching accurate results in terms of free energy estimates.
Here, we first derive a simple metadynamics reweighting scheme starting from the approach proposed by Tiwary et al. (16) and taking the limit of short trajectories. We then analyze its properties compared to Tiwary’s approach. Then, we relate the new scheme to the one from Bonomi et al. (14) After that, we numerically evaluate the accuracy of the different reweighting schemes to reproduce the underlying free energy landscape as a function of the simulation time and compare it to the one obtained from the bias potential. We do this first using a simple unidimensional system and then the alanine dipeptide, which represents a standard benchmark for these calculations. The latter case is analyzed both in the context of standard and well-tempered metadynamics.

2. Results and Discussion

ARTICLE SECTIONS
Jump To

2.1. Reweighting Models

In a standard metadynamics simulation the time-dependent bias potential V(r,t) at time t and at point r in coordinate space is defined as a function of the past trajectory of the system. In particular it is a sum of Gaussian hills defined in the d-dimensional space of collective variables s with height W and widths σ, deposited with a period Δt:
(1)
where are the Gaussian hills. The system evolves according to the total internal energy U(r) + V(s(r),t), where U(r) is the unbiased internal energy of the system.
The conformations of the system are sampled along the trajectory using the time-dependent biased energy function. Eventually, however, we want to obtain the properties (i.e., averages of observables) of the unbiased system. In the unbiased system, the statistical weight of each sampled conformation is the same. In metadynamics simulations, each sampled conformation is the result of an increasing bias potential, so its statistical weight w(r,t) is also a function of the bias potential. A right choice of w(r,t) insures that unbiased averages of an arbitrary operator A can be calculated via
(2)
with ⟨···⟩ 0 denoting the unbiased average, and ⟨···⟩b denoting the average over the biased trajectory.
For a system biased with an external potential V(s(r),t) the instantaneous probability density function (pdf) is given as
(3)
The corresponding unbiased system has the pdf:
(4)
To recover unbiased averages from a biased simulation, the trajectory averages have to be weighted according to eq 2. The statistical weight can be calculated as follows: (14,16)
(5)
(6)
Since the bias potential only depends on r indirectly through s, the same is the case for w, and the expression above can be further modified: (14,16)
(7)
(8)
(9)
where
(10)
is the unbiased pdf in the low-dimensional CV space Ω and F(s), the unbiased free energy.
Metadynamics provides a way to estimate the unbiased free energy F(s) from the bias potential (12,13,15,22)
(11)
where γ is the so-called bias factor of well-tempered metadynamics (in standard metadynamics γ → ) and c(t) is a time-dependent offset which can be chosen so that F(s) becomes asymptotically time-independent. (16)
Tiwary and Parrinello (16) propose to replace the instantaneous estimate of F(s) from eq 11 into eq 9, which leads to Tiwary’s weights. (16,19) In the case of standard metadynamics, this results in
(12)
where the integral in ds is extended to the sampled regions of the collective variable space. Thus, the weight according to this scheme is proportional to the exponential of the bias potential normalized by the average of the exponential of the bias potential over the space of the collective variables.
Alternatively, we could replace into eq 9 the more accurate estimate of F(s) from the end of the simulation, resulting in the following expression (for standard metadynamics):
(13)
Here tf is the final time of the simulation. In the later stages of the trajectory, where V(s,t) → V(s,tf), eq 13 converges to the Tiwary’s weights of eq 12. On the other hand, considering for simplicity the case of standard metadynamics, in the early part of the trajectory for ttf, V(s,tf) dominates in the term at the numerator, so the integrals at the numerator and denominator of eq 13 are approximately equal, and the expression simplifies to exponential weights:
(14)
Thus, eq 14 and eq 12 represent the two limiting cases of eq 13 for short and long trajectories, respectively. It is worth noting that, in our tests, the reconstructed FES obtained with the reweighting scheme of eq 13, although getting close to the accuracy of the exponential reweighting and to the Tiwary reweighting in short and long simulations, respectively, showed always slightly lower accuracies than the other two, so we will include in our analysis only the Tiwary eq 12 and a modified version of the exponential reweighting as detailed below.
The plain exponential reweighting as in eq 14 may provide enhanced accuracy on short standard metadynamics trajectories, however it has the disadvantage that the weight of conformations sampled early along the trajectory will be vanishingly small with respect to those sampled at a later stage, due to the constantly increasing bias potential in the exponential. To obviate this problem, we noted that the same metadynamics trajectory would be generated if, beside regularly adding Gaussians gt(s) to the bias potential, a constant term would be subtracted from it, uniformly over the space of s. In this way, we will have a modified bias potential
(15)
where is the average of the total bias potential obtained summing up the volume of all the Gaussians until time t. The modified bias has the property that at any time t. Now, replacing the bias V in eq 14, with V′ from eq 15, leads to a balanced exponential reweighting scheme:
(16)
This scheme can be implemented very simply without any modification of the available software for metadynamics simulations. Comparing the balanced exponential reweighting eq 16 and Tiwary’s reweighting eq 12 shows that they differ in the way the leading exponential term is normalized. In Tiwary’s method, it is normalized with respect to ⟨eβV(s,t)s, while in the balanced exponential, it is normalized by eβ⟨V(s,t)⟩s. The average of the exponential in Tiwary’s method, although analytically more accurate, converges more slowly than the exponential of the average bias because it is very sensitive to values of V at the upper tail of its distribution. For the same reason, Tiwary’s method may result in a larger run-to-run variability due to the different ways the distribution of V may be progressively sampled in different runs.
In what follows, we will show how the balanced exponential reweighting scheme in eq 16 relates to the method by Bonomi et al., (14) and then, in the next sections, we will test it. In the scheme proposed by Bonomi et al. (14) the histogram of the collective variables is accumulated in between two Gaussian deposition steps, and then it is reweighted according to
(17)
where is proportional to the Gaussian being added to the bias potential and is proportional to its biased average over the CV space. In the original publication, (14) the biased average is obtained directly from the biased histogram. In the program distributed with PLUMED 1.3, which implements this algorithm, another method is also available (with the option -rewtype 0, whereas the original method corresponds to -rewtype 1) the biased average is obtained by averaging over the estimated biased distribution Pb(s,t) = eβV (s,t)/(γ–1). In our hands the second method provides more accurate estimates of the free energy in most of the circumstances, so we will only show those results. From the biased histogram, the unbiased histogram is obtained as P0(s) = eβV (s,tf)Pb(s,tf). This method is based on the construction and time evolution of histograms, rather than the assignment of weights to sampled conformations. To cast it into weights, let us define nh(s) the partial histogram obtained from the data collected between the deposition of the (h – 1)th and hth Gaussians. After deposition of the hth Gaussian, the partial histogram is added to overall histogram and reweighted according to eq 17. Thus, the weight of the partial histogram on the overall histogram becomes nh(s)e–β((s,t)+ċ(t))Δt. At the next Gaussian deposition, a similar process occurs, and the weight of our initial partial histogram becomes
(18)
The same process is repeated for all subsequent Gaussians added to the bias potential. The weight of the partial histogram on the final overall unbiased histogram is then obtained as
(19)
Here we used the fact that V(s, th) = = . Equation 19 tells us that every sample collected between the (h – 1)th and the hth Gaussian deposition has a weight:
(20)
Next, we note that, in the exponential expression for the weight in eq 20, there are M terms, where M is the total number of Gaussians added to the bias at the final time t = tf. The first h – 1 terms, all included in the expression V(s,th–1), come from the Gaussians deposited until time th–1 and are dependent on s; the rest are uniform offsets coming from Gaussians deposited at later stages. The expression in the exponential of eq 20 is defined apart from an additive constant C. If we choose we obtain
(21)
where ⟨−⟩b indicates an average over the biased ensemble. We note that eq 21 recalls eq 16 with the difference that in eq 16 the averages of the Gaussians are simple averages while in eq 21 they are made in the biased ensemble. In the case of standard metadynamics, where the bias factor γ → , the averages in the biased distribution eβV(s,t)/(γ–1) become simple averages in s space and the two methods should provide approximately similar results, apart from the fact that in Bonomi’s methods all evaluations of the bias potential are done on binned data, rather than the unbinned data used in eq 16. This latter difference turns out to have important consequences in the case of estimating averages of observables that were not directly biased during the metadynamics simulations, as will be shown later. In any case, for the comparison, we will use the output of the program from ref (14) made available via PLUMED 1.3 and not eq 20.
Another scheme that is often used is the one from Branduardi et al. (15) which assumes that the weights only depend on the final value at time tf of the bias potential and are not a function of time.
(22)
This is suitable for well-tempered metadynamics, where the time-dependent changes to the bias potential, after a transient, become exponentially small, and so the bias becomes effectively constant.

2.2. Unidimensional Model System

The effects of the reweighting models specifically indicated for standard metadynamics, i.e., the balanced exponential (eq 16) and Tiwary’s (eq 12) reweighting scheme, and presented above have been tested on a simple toy model: a particle in the double well potential U(x) = (x2 – 1)2. The simplicity of the model allows us to associate precisely the effects of each parameter of the metadynamics to the dynamics of the particle. For the sake of testing the reweighting models, we do not even need to run a real molecular dynamics. Instead we built a simple Monte Carlo scheme where the particle can move randomly in the interval [x – Δx: x + Δx] where x is its present position and Δx is fixed to 0.05. Moves are accepted according to a Metropolis scheme where the energy of the system is U(x) + V(x,t), V(x,t) being the metadynamics bias.
(23)
Here Δt is the deposition period, v is the volume, and σ is the width of the Gaussian hills (the height of the Gaussian is then ). We explored a wide spectrum of values of the three metadynamics parameters, while we fixed the length of the simulations to 2 × 107 steps. The temperature of the system used in the Metropolis scheme was set to 1/10th of the height of the free energy barrier separating the two minima, which insures that unbiased simulations would not converge quickly. A total of 72 independent trajectories were collected for each set of parameters. The position of the particle was saved every 100 steps unless reported differently and, later, discretized into bins of size 0.005.
The true free energy landscape of the single particle in the double-well potential corresponds trivially to the potential F(x) = U(x) itself. We used the simulation data to estimate it either using the negative bias potential
(24)
or the position distribution reweighted according the two reweighting schemes suitable for standard metadynamics presented above (eqs 12 and 16):
(25)
Here the δ functions are to be considered as characteristic functions of the corresponding discretized position bin. We note that for standard metadynamics simulations replacing the final negative bias potential with its time average over the last part of the simulation, after the relevant CV space has been sampled, can provide a more accurate approximation of the free energy landscape. (23,24) The averaging procedure would produce its effects mostly on the free energy estimates at the late stages of the simulation when enough statistics have been accumulated. In addition, the same averaged potential could be used in the reweighting algorithms, instead of the instantaneous potential, producing a similar improvement of the derived free energy estimates. Since we are interested in the relative behavior of the free energy estimates by the various schemes and in particular in the speed of convergence in the early stages of the simulation, we have not extended our analysis to the case of the time averaged bias potential.
Examples of estimated free energy landscapes for one particular simulation are reported in (Figure 1a). We characterized the deviations of the estimated FESs from the true one using the error or root-mean-squared deviation , where nx is the number of discretized position bins in the sum. The sum is limited to the values of x sampled during the simulation and to the interval [-1.5,1.5] to exclude the external parts of the landscape which may reach very high free energy values and provide dominant contributions to the integral. The FESs obtained with the various methods are aligned by subtracting the average value over the points in the sum to obtain the RMSDFES. The RMSD values are averaged across all runs with the same metadynamics parameters. The error of the free energy estimates measured at the end of the simulations is high at low filling/flooding rates vt, where the length of the simulations is not sufficient to properly sample the free energy landscape, as there is no time to cross the free energy barrier. For sufficiently large flooding rates, instead, the error becomes approximately proportional to as already reported. (25) We then fixed the metadynamics parameters to values providing a sufficiently fast flooding rate and monitored the average error as a function of simulation time (Figure 1b). In this case the balanced-exponential reweighted FES and, to a lower extent, the negative bias potential converge significantly faster to the true FES than the Tiwary reweighted FES. The balanced exponential FESs show run-to-run variability lower than Tiwary’s method. Both balanced exponential and Tiwary’s method show that the tends to a constant toward the end of the simulations (Figure 1b inset), as expected. (14) The accuracy of the negative bias reaches a plateau, related to the continuous addition of new hills of fixed height to the bias potential, as per standard metadynamics, which produces finite fluctuations around the true FES. The fluctuations of the balanced-exponential and the Tiwary reweighted FES around the correct FES decrease with time due to the time dependence of the size of the weights (see below). The behavior of the fluctuations of the two reweighted FESs can also be observed directly in the Figure 1a.

Figure 1

Figure 1. (a) FES obtained along one run after 106, 2 × 106, 3 × 106, 4 × 106 and 2 × 107 steps using balanced exponential reweighting, Tiwary reweighting and negative bias potential (red, blue and green points, respectively). The true potential (purple) is superimposed for comparison. (b–e) Time series of (b) the RMSDFES (in inset ), (c) the absolute value of the estimated free energy difference between the two minima of the double well potential, (d) the absolute value of the deviation of the estimated height of the free energy barrier from the true value of 1, and (e) the RMSDV. Same color scheme as in part a. The solid lines represent the average values of the quantities across the 72 independent runs. Shaded bands indicate the standard deviations. The positions of free energy minima and maxima are computed by fitting parabolas to the data points available in the intervals [−1.2,0.8], [−0.2,0.2], and [0.8,1.2] and taking the ordinate of the vertexes. (f–g)Time series of the position of the particle along one run where balanced exponential (f) and Tiwary (g) weights are reported according to a color scale.

Inspection of the estimated FESs from one of the simulations (Figure 1a) shows that in the early stages of the simulation, when the negative bias has not converged yet in the barrier region because of limited sampling, the balanced-exponential reweighted FES provides already an estimate close to the true FES. In other words, the reweighting helps to correct for the limited sampling in the high free energy regions. The time series of the free energy difference (|ΔF|) between the two minima of the double well potential averaged across the independent runs (Figure 1c) shows that the bias potential and the balanced-exponential reweighting reach low values faster than the method proposed by Tiwary and Parrinello, although the fluctuations in this case are larger than those for the RMSD. The balanced exponential and Tiwary’s method eventually converge to similar values smaller than the negative bias and with smaller fluctuations. The time series of the height of the free energy barrier ΔF (or, better, its deviations from the expected value of 1, Figure 1d), which measures how well the methods deal with high free energy regions like transition states that are sampled rarely, shows, on average, that the balanced-exponential reweighted FES manages to converge several times faster than the Tiwary reweighted one or the bias potential, and it does so consistently across the runs (small standard deviation). Also in this case, eventually, on longer simulation time scales, Tiwary’s method reaches a similar accuracy.
By analyzing the behavior of the weights along the simulations, we observed that in the balanced exponential reweighting method (Figure 1f), the weights remain relatively constant although early points show slightly lower weights. This is in contrast with what would happen using a plain exponential reweighting scheme where the weights would increase exponentially along the trajectory resulting in early points being “forgotten” due to low weights compared to later points. The weights in the Tiwary method (Figure 1g) become quickly constant in each region of coordinate space that is visited, although early points show slightly larger weights than subsequent ones. The constancy of weights for similar coordinates at the beginning and the end of the simulation insures that all points observed along the simulations contribute to the averages and eventually lead to a better convergence to the underlying FES. The differences in the weights of the initial part of trajectory in Tiwary’s and balanced-exponential scheme may be at the origin of the different convergence speed of the two methods. A careful analysis of the weights across the runs and how each data point contributes to the averages will be provided later for the case of the alanine dipeptide.
In the simple model analyzed here, we can compare the results directly with true FES, but in the normal case, the underlying FES is not known. However, we can compare the reconstructed FES with the negative bias potential, which represents a well-defined reference. The comparison shows that, in the early stage of the simulation, the deviation () between negative bias and balanced-exponential reweighted FES is smaller than the one between the negative bias and the Tiwary’s reweighted FES (Figure 1e). Then, the two lines converge to similar values in analogy to what is seen for RMSDFES. Compared to parts b–e of Figure 1, it seems that the convergence of RMSDV of the two reweighting techniques to within a standard deviation occurs at a time larger than required for the balanced-exponential FES to converge to the true FES. The time of convergence of the RMSDV of the two reweighting schemes may then provide an indication about the state of convergence of the simulation.

2.3. Alanine Dipeptide with Standard Metadynamics

After having tested our findings on a simple 1-dimensional system, we applied the same analysis to a more complex system, the alanine dipeptide in vacuum, which has been used in the literature as a standard benchmark. (14−16) In this case, the biasing variables are the Ramachandran angles φ and ψ characterizing the conformation of the peptide. A 75 ns-long well-tempered metadynamics simulation of the peptide has been used to retrieve an accurate reference free energy landscape as a function of the Ramachandran angles (Figure 2a), using 1.2 kJ/mol, 0.35 rad, and 6 for the hills’ height, width, and bias factor, respectively. Hills were deposited every 1.0 ps. The simulation were carried out using the program GROMACS 2018.1 (26) and the PLUMED plugin. (19) The Amber03 force field (27) was used resulting in a free energy landscape showing a well-defined set of features which include the minimum C7eq (φ = −1.382, ψ = 1.005), the local minimum Cax (φ = 1.257, ψ = −0.880), as well as the connecting transition state ‡ (φ = 1.885, ψ = −2.136).

Figure 2

Figure 2. (a) Reference FES of the alanine dipeptide as a function of the CV φ and ψ. Contour levels are plotted every 1 kbT. The black contours indicate the region within 10 kbT from the global minimum, which is used for the RMSDFES calculations. The red triangle, square and circle indicate the position of the minima C7eq and Cax and the transition state ‡, respectively. The color shades are guides for the eye. (b–f) Difference between the reference FES and those estimated using the balanced exponential, Tiwary’s, Branduardi’s, and Bonomi’s reweighting, and the negative bias, respectively, after 1.4 ns in one of the runs. Contour levels are plotted every 0.25 kbT, and the range of color shades is 4 times smaller than that in part a.

To test the reweighting schemes, 96 standard metadynamics simulation runs were performed using 1.004 kJ/mol, 0.25 rad, and 2 ps as Gaussian height, width, and deposition period, respectively, similar to ref (16). Standard metadynamics was obtained by setting the bias factor to a large value (γ = 60000). The reweighting factor was updated after every hill deposition. The values of the collective variables were saved every 10 MD steps. In all these simulations, we estimated the FES using the balanced exponential, the Tiwary reweighting schemes, and the negative bias potential as a function of the simulation time and compared them to the reference FES. Although not specifically indicated for standard metadynamics, we also included Bonomi’s (14) and Branduardi’s (15) methods to the comparison. The reference FES and the deviations of the FES estimates for the various reweighting schemes along one run after 1.4 ns simulation are reported in Figure 2. To quantify the differences, we measured the root mean squared deviation of the reweighted FESs (after aligning them as for the monodimensional case) from the accurate reference (, where the sum is extended to the regions within 10kbT from the global minimum and nφ,ψ is the number of bins used in the sum, see Figure 3a), as well as the deviation |ΔΔFC7eqCax| of the estimated free energy difference between minima from the reference (Figure 3b), and the deviation |ΔΔFC7eq–‡| of the estimated height of the free energy barrier (Figure 3c).

Figure 3

Figure 3. (a–c) Time series of the RMSDFES and the deviation from the reference ΔFC7eqCax and ΔFCax–‡, respectively, for the balanced exponential (red), Tiwary’s (blue), Branduardi’s (purple), and Bonomi’s (orange) reweighting, and negative bias (green) estimate of the FES. The data are averaged over 96 runs. Shaded bands indicate standard deviations. The black line at 1.4 ns indicates the time point where the FESs in Figure 2 have been extracted. (a, inset) shows to approximately reach a plateau for balanced exponential, Tiwary’s, and Bonomi’s reweighting schemes. (d) RMSDV between the different reweighting schemes and negative bias (same color scheme as above). (e, f) Time series of ψ when φ = −1.88 rad along one selected run. Balanced exponential (e) and Tiwary (f) weights are reported according to the color scale. ⟨w⟩ is the average of the weights along the run.

The data show that, in the early stages of the simulation, the balanced-exponential and, to a slightly lower extent, Bonomi’s reweighted FES provide smaller differences from the accurate FES than the other methods (Figure 3a), as already observed in the 1-dimensional case. Also the run-to-run variability of balanced exponential and Bonomi’s method is smaller than the other methods as shown by the width of the error bands in Figure 3a. Only in the later stages does the Tiwary reweighted FES become, on average, as accurate as the other two methods mentioned above, although this method is affected by a very large run-to-run variability, with a standard deviation of the RMSDFES that is 3 times larger than the balanced exponential. Closer inspection shows that in the runs where the local minimum Cax is reached early in the simulation, Tiwary’s reweighting does not converge on the reference FES due to wrong weighting of the conformations of the local minimum (Figures S1 and S2 in the Supporting Information).
Also, the size of the weights along the simulation trajectory (Figure 3ef) shows a strong analogy to the 1-dimensional case: in the balanced exponential reweighting, the weights are slightly damped in the early part of the trajectory, while in the Tiwary reweighting early points are slightly heavier. To evaluate more carefully the effect of the reweighting schemes, we measured the Kish’s effective sample size () (28) in the case of the balanced exponential and Tiwary’s method (Figure 4). The effective sample size for the balanced exponential reweighting shows a slow down in its growth after approximately 200 ps, and then it starts to grow again, while Tiwary’s sample size continues to grow steadily although with a very small slope. In the interesting regime around 1 ns of simulation (sample size 5 × 104), where the exponential reweight outperforms Tiwary’s reweight in terms of matching the reference free energy, the effective sample size of Tiwary’s and balanced exponential reweighting are approximately 0.19 and 0.06 of the total sample size, respectively, so less than an order of magnitude difference. This means that the lower size of the sample in balanced exponential reweighting is more than compensated for by the quality of the weights in terms of matching the reference FES.

Figure 4

Figure 4. Kish’s effective sample size of the Tiwary’s (purple) and balanced exponential (green) reweighting method as a function of the total sample size along one standard metadynamics trajectory of the alanine dipeptide (data saved every 0.02 ps). The blue line represents the size of the ideally uncorrelated sample.

The negative bias and Branduardi’s FES are changed by a finite and constant amount every time a new hill is added to the bias as per standard metadynamics, so they show finite fluctuations around the reference FES which do not decrease in size over time. Balanced exponential, Tiwary’s and Bonomi’s methods, on the other hand, also thanks to a more uniform distribution of weights along the simulations as mentioned above, allow for a continuous improvement of the FES estimate similar to the one expected for well-tempered metadynamics where (Figure 3a, inset).
Similar to the 1-D case, the difference with the negative bias is smaller for the balanced exponential FES than Tiwary’s method in the early stages of the simulation (Figure 3d). Bonomi’s method follows balanced exponential closely. At the later stages Tiwary’s, Bonomi’s, and balanced exponential reweighting converge on a similar RMSDV. Also, in this case, the time of convergence of Tiwary’s and balanced exponential provides and indication of the degree of convergence of the simulation. In this regime, Branduardi’s methods shows a much smaller distance to the negative bias than the other methods. Branduardi’s method indeed is supposed to asymptotically match the negative bias.
The data provided show that Branduardi’s reweighting method cannot be productively applied straight away to standard metadynamics data, since it assumes that the bias potential does not change significantly anymore after a transient, which is only the case in well-tempered metadynamics. It could be applied to the time averaged bias potential, where the fluctuations due to hills addition decrease with time. (30) The data show that although RMSDFES remains higher than the other methods, the lower parts of the FES are reproduced with relative accuracy (low error on the free energy difference between minima) while the high free energy regions are not (large error on the height of the free energy barrier). On the other hand, Bonomi’s method, which is also not specifically indicated for standard metadynamics, provides good results also in this context, with accuracies similar to those of the balanced exponential method, confirming the similarity of the methods in the case of standard metadynamics highlighted in the section Reweighting Models.
The data presented above show the variability of the FES estimates across independent runs, thus providing a stringent measure of the accuracy of the different approaches. It is however also important to be able and provide an estimate of the accuracy of a FES directly from a single run using block analysis. (29) We have performed a block analysis of a randomly selected standard metadynamics simulation of the alanine dipeptide. The results (Figure 5) show that the errors on the reweighted FES obtained using the balanced exponential reweighting are similar, although only slightly higher than those obtained by Tiwary’s reweighting method, in the regime of simulation lengths we are interested in.

Figure 5

Figure 5. Block analysis of the average error on the determination of the free energy landscape of the alanine dipeptide using standard metadynamics simulations and the Tiwary’s (green) or the balanced exponential (purple) reweighting scheme. The computation is performed on a randomly selected 2 ns-long trajectory. The error is averaged over the regions of the FES where the free energy is within 10kbT of the global minimum.

2.4. Alanine Dipeptide with Standard Metadynamics on a Single Collective Variable

Reweighting specifically allows for evaluating averages of variables that were not used in biasing the simulation. To evaluate the accuracy of this procedure with the various reweighting schemes, we conducted standard metadynamics simulations as done in the previous section but biasing only the φ dihedral angle (and decreasing by a factor of 10 the height of the Gaussians, which helps to approximately match the simulation time required to fill up the minima in this case). Then, we evaluated the 2-dimensional FES as a function of both φ and ψ using the various reweighting schemes and compared the results to the reference free energy as done before. This represents a sort of worst-case scenario where the secondary variable ψ represents an almost independent degree of freedom with respect to φ.
The deviations of the reweighted FESs from the reference in a typical run are reported in Figure 6 (Figure S3 in the Supporting Information, for a run sampling the local minimum Cax early in the trajectory). In this case, a major problem is the limited sampling in the ψ direction, which, as mentioned, has low correlation with the φ direction used in the bias. In particular, for a given φ “slice” of the FES corresponding to a minimum of the free energy, the values of ψ with high free energy will be sampled rarely, because of the bias pushing the system away from there based only on the value of φ. This problem appears in the balanced exponential reweighting scheme, where outliers are present (Figure 6b) with very low weights (i.e., very high free energy) located at the border of the sampled regions and corresponding to regions visited early in the simulation and never visited again. These outliers contribute substantially to increase the RMSDFES although in reality the rest of the reweighted FES is close to the reference. To limit the influence of this phenomenon, we measured the RMSDFES of the various reweighting schemes limiting the sum of the differences to the bins visited in the last 1 ns of simulation, thus excluding the few bins with unusually high free energy due to limited sampling.

Figure 6

Figure 6. Same as Figure 2 but for the standard metadynamics runs using only φ as biasing collective variable.

The RMSDFES (Figure 7a), in this case, shows that balanced exponential and Tiwary’s methods are, on average, more accurate than the others, both in the early and the later stages of the simulations followed by Branduardi’s method. However, Tiwary’s and Branduardi’s method are both affected by a large run-to-run variability: the standard deviation of the RMSDFES is about twice as big as the balanced exponential (width of the error bands in Figure 7a). In the case of Tiwary’s method, this is again due to a fraction of the runs entering the local minimum Cax early in the trajectory (Figure S3 in Supporting Information). Also the estimate of the free energy difference between minima in the case of the balanced exponential (Figure 7b) converges faster to more accurate values than Tiwary’s and Bonomi’s methods and shows smaller run-to-run fluctuations. The data on the estimate of the height of the free energy barrier are to be taken with care, as the TS lies close to the border of the sampled region (Figure 6) and suffers from limited sampling. In any case, the balanced exponential shows on average to estimate the barrier accurately followed by Tiwary’s and Branduardi’s methods. Here we observed that, in Bonomi’s method, using the averaging of the hills over the biased distribution (rewtype = 0) does not provide very good estimates of the free energy barrier height. Averaging by means of the histogram (rewtype = 1) showed better results, similar to Branduardi’s method (data not shown).

Figure 7

Figure 7. Same as Figure 3 but for the standard metadynamics runs using only φ as biasing collective variable.

Although Bonomi’s and balanced exponential reweighting should provide very similar weights in standard metadynamics simulations, the two methods differ in how the averages of observables not biased in the metadynamics are treated. Bonomi’s method relies on the time evolution of the histogram, while the balanced exponential creates a histogram of weighted data. This difference may be at the origin of the unexpected discrepancies observed in the data shown above. We also note that, as already reported in ref (14), the accuracy of Bonomi’s method may be influenced by the bin size and sampling rate, which in the present case (0.126 rad and 20 fs, respectively) do not match the optimal values used in that publication.

2.5. Alanine Dipeptide with Well-Tempered Metadynamics

As a final test, we performed the same analysis on a set of 96 well-tempered (WT) metadynamics runs of the alanine dipeptide in vacuum. In this case, the same metadynamics parameters as for the construction of the reference FES have been used (see above). We evaluated the same quantities identified in the previous sections to determine the accuracy of the reconstructed FESs (Figure 8). The data show a slightly different picture with respect to standard metadynamics. The estimated FESs show a similar accuracy in reproducing the correct free energy difference between minima (Figure 9b), as well as the height of the free energy barrier in the early stages of the simulations (Figure 9c), with the exception of Bonomi’s method, which takes slightly longer to converge. In terms of RMSDFES (which is focused on the regions around the minima), the balanced exponential and Tiwary’s estimates are only slightly more accurate than the other methods, on average (Figure 9a), although the first one shows lower run-to-run variability than the latter (i.e., standard deviation of RMSDFES is almost half) in analogy to the previously analyzed cases. The Branduardi estimate, after a transient, outperforms rescaled negative bias estimates (as already observed in ref (15)). Bonomi’s reweighting requires a slightly longer transient to reach accuracies similar to the other methods but in the later stages of the simulations there is substantial overlap among the various methods, confirming what reported in ref (21) for longer time scales. The fact that the height of the Gaussians is reduced during the course of the simulations helps to reduce the weight differences between early and late trajectory snapshots (Figure 9e,f) and insures that the for all FES estimates (Figure 9a inset and Figure S4 in the Supporting Information). One other aspect that is similar to both standard and WT-metadynamics is that, in the early stage of the simulations, the estimated FES from balanced exponential reweighting is closer to the estimate obtained from the negative bias than Tiwary’s estimate (Figures 9d and 3d), and the convergence of the two values can be used as an indication of the convergence of the simulation. As a final note, we would like to point out that the estimates obtained using the balanced exponential, Bonomi’s and Tiwary’s method on standard metadynamics runs reach (after 8 ns) accuracies comparable to those obtained using the rescaled negative bias from WT-metadynamics runs (Figure 9a–d and 3a–d).

Figure 8

Figure 8. Same as Figure 2 but for the well-tempered metadynamics runs.

Figure 9

Figure 9. Same as Figure 3 but for the well-tempered metadynamics runs. In the inset of panel a, the reaches a plateau, as expected, although this is more evident when using a larger free energy cutoff in the calculation of the RMSDFES as shown in Figure S4 in the Supporting Information.

3. Conclusions

ARTICLE SECTIONS
Jump To

In this work, we analyze the behavior of several reweighting schemes for metadynamics simulations. We specifically focus on the early stages of the simulations and the speed of convergence of the methods. We show that, although the free energy surface obtained by reweighted data can be more accurate than the negative bias potential, none of the previously published reweighting schemes tested here provides optimal results in all the situations we have considered. Tiwary’s method (16) works well in well-tempered metadynamics but converges slowly in the case of standard metadynamics and in all cases shows large run-to-run variability due to overweighting early trajectory snapshots. Bonomi’s method (14) converges quickly and accurately in the case of standard metadynamics but does not do the same in well-tempered metadynamics and, when dealing with not-biased observables, does not provide results as accurate as other methods.
We have then introduced a balanced exponential reweighting scheme that can be derived in the limit of short trajectories from a modified Tiwary’s reweighting scheme. We have highlighted the reasons why the new scheme converges faster to the underlying free energy than Tiwary’s method and with smaller run-to-run fluctuations, and we have shown that it is closely related to Bonomi’s approach in the case of standard metadynamics. The numerical tests show that the free energy surface reconstructed using the balanced exponential scheme, on average, converges in all the tested situations to the underlying free energy surface in the same time or faster than the best of the other analyzed reweighting schemes and with low run-to-run variability.
In addition, we show that the accuracy provided by the balanced exponential and Tiwary’s and Bonomi’s reweighted FES in the later stages of standard metadynamics simulations may be comparable with the one obtained by well-tempered metadynamics, where the effect of reweighting on the accuracy of the FES is minor. Proper sampling with well-tempered metadynamics is dependent on the choice of the bias factor, which involves an estimate of the height of the free energy barriers of the investigated system. Since this estimate may not always be possible due to system size and complexity, standard metadynamics with reweighted FES may provide a valid alternative. In those cases, it is often not feasible to continue simulations past the point of convergence, so reweighting has to be accurate even at early stages of the simulation to make efficient use of the limited amount of data. The proposed balanced exponential reweighting scheme may be able to address those situations, especially in cases where the biased collective variables capture the major degrees of freedom of the system under scrutiny.

Supporting Information

ARTICLE SECTIONS
Jump To

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.9b00867.

  • Figures S1–S4 (PDF)

Terms & Conditions

Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system: http://pubs.acs.org/page/copyright/permissions.html.

Author Information

ARTICLE SECTIONS
Jump To

  • Corresponding Author
  • Author
    • Timo M. Schäfer - Institut für Physik, Johannes Gutenberg-Universität Mainz, Mainz, GermanyGraduate School Materials Science in Mainz, Mainz, Germany
  • Notes
    The authors declare no competing financial interest.

Acknowledgments

ARTICLE SECTIONS
Jump To

We thank Giovanni Bussi for critical reading of the manuscript and Friederike Schmid and Marialore Sulpizi for fruitful discussions. T.M.S. gratefully acknowledges financial support from the Graduate School Materials Science in Mainz. G.S. gratefully acknowledges financial support from the Max-Planck Graduate Center with the University of Mainz. We gratefully acknowledge support with computing time from the HPC facility Hazelhen at the High Performance Computing Center Stuttgart and the HPC facility Mogon at the university of Mainz. This work was supported by the German Science Foundation within the SFB 1066 Project Q1.

References

ARTICLE SECTIONS
Jump To

This article references 30 other publications.

  1. 1
    Bernardi, R. C.; Melo, M. C.; Schulten, K. Enhanced sampling techniques in molecular dynamics simulations of biological systems. Biochim. Biophys. Acta, Gen. Subj. 2015, 1850, 872877,  DOI: 10.1016/j.bbagen.2014.10.019
  2. 2
    Dror, R. O.; Dirks, R. M.; Grossman, J.; Xu, H.; Shaw, D. E. Biomolecular Simulation: A Computational Microscope for Molecular Biology. Annu. Rev. Biophys. 2012, 41, 429452,  DOI: 10.1146/annurev-biophys-042910-155245
  3. 3
    Abrams, C.; Bussi, G. Enhanced sampling in molecular dynamics using metadynamics, replica-exchange, and temperature-acceleration. Entropy 2014, 16, 163199,  DOI: 10.3390/e16010163
  4. 4
    Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23, 187199,  DOI: 10.1016/0021-9991(77)90121-8
  5. 5
    Hansmann, U. H. Parallel tempering algorithm for conformational studies of biological molecules. Chem. Phys. Lett. 1997, 281, 140150,  DOI: 10.1016/S0009-2614(97)01198-6
  6. 6
    Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 1999, 314, 141151,  DOI: 10.1016/S0009-2614(99)01123-9
  7. 7
    Sugita, Y.; Okamoto, Y. Replica-exchange multicanonical algorithm and multicanonical replica-exchange method for simulating systems with rough energy landscape. Chem. Phys. Lett. 2000, 329, 261270,  DOI: 10.1016/S0009-2614(00)00999-4
  8. 8
    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, 291318,  DOI: 10.1146/annurev.physchem.53.082301.113146
  9. 9
    Darve, E.; Rodríguez-Gómez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations. J. Chem. Phys. 2008, 128, 144120,  DOI: 10.1063/1.2829861
  10. 10
    Maragliano, L.; Vanden-Eijnden, E. A temperature accelerated method for sampling free energy and determining reaction pathways in rare events simulations. Chem. Phys. Lett. 2006, 426, 168175,  DOI: 10.1016/j.cplett.2006.05.062
  11. 11
    Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 1256212566,  DOI: 10.1073/pnas.202427399
  12. 12
    Bussi, G.; Laio, A.; Parrinello, M. Equilibrium free energies from nonequilibrium metadynamics. Phys. Rev. Lett. 2006, 96, 090601  DOI: 10.1103/PhysRevLett.96.090601
  13. 13
    Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100, 020603,  DOI: 10.1103/PhysRevLett.100.020603
  14. 14
    Bonomi, M.; Barducci, A.; Parrinello, M. Reconstructing the equilibrium Boltzmann distribution from well-tempered metadynamics. J. Comput. Chem. 2009, 30, 16151621,  DOI: 10.1002/jcc.21305
  15. 15
    Branduardi, D.; Bussi, G.; Parrinello, M. Metadynamics with adaptive Gaussians. J. Chem. Theory Comput. 2012, 8, 2247,  DOI: 10.1021/ct3002464
  16. 16
    Tiwary, P.; Parrinello, M. A Time-Independent Free Energy Estimator for Metadynamics. J. Phys. Chem. B 2015, 119, 736742,  DOI: 10.1021/jp504920s
  17. 17
    Tiana, G. Estimation of microscopic averages from metadynamics. Eur. Phys. J. B 2008, 63, 235,  DOI: 10.1140/epjb/e2008-00232-8
  18. 18
    Smiatek, J.; Heuer, A. Calculation of free energy landscapes: A histogram reweighted metadynamics approach. J. Comput. Chem. 2011, 32, 20842096,  DOI: 10.1002/jcc.21790
  19. 19
    Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185, 604613,  DOI: 10.1016/j.cpc.2013.09.018
  20. 20
    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 dynamics. Comput. Phys. Commun. 2009, 180, 1961,  DOI: 10.1016/j.cpc.2009.05.011
  21. 21
    Gimondi, I.; Tribello, G. A.; Salvalaglio, M. Building maps in collective variable space. J. Chem. Phys. 2018, 149, 104104,  DOI: 10.1063/1.5027528
  22. 22
    Dama, J. F.; Parrinello, M.; Voth, G. A. Well-tempered metadynamics converges asymptotically. Phys. Rev. Lett. 2014, 112, 240602,  DOI: 10.1103/PhysRevLett.112.240602
  23. 23
    Laio, A.; Gervasio, F. L. Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys. 2008, 71, 126601,  DOI: 10.1088/0034-4885/71/12/126601
  24. 24
    Marinelli, F.; Pietrucci, F.; Laio, A.; Piana, S. A Kinetic Model of Trp-Cage Folding from Multiple Biased Molecular Dynamics Simulations. PLoS Comput. Biol. 2009, 5, e1000452  DOI: 10.1371/journal.pcbi.1000452
  25. 25
    Laio, A.; Rodriguez-Fortea, A.; Gervasio, F. L.; Ceccarelli, M.; Parrinello, M. Assessing the Accuracy of Metadynamics †. J. Phys. Chem. B 2005, 109, 67146721,  DOI: 10.1021/jp045424k
  26. 26
    Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A message-passing parallel molecular dynamics implementation. Comput. Phys. Commun. 1995, 91, 4356,  DOI: 10.1016/0010-4655(95)00042-E
  27. 27
    Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber biomolecular simulation programs. J. Comput. Chem. 2005, 26, 16681688,  DOI: 10.1002/jcc.20290
  28. 28
    Kish, L. Survey Sampling; John Wiley & Sons, Ltd.; New York, 1965.
  29. 29
    Frenkel, D.; Smit, B. Understanding Molecular Simulations, 2nd ed.; Computational Science; Academic Press: 2002; Vol. 1.
  30. 30
    Bussi, G. Personal communication.

Cited By

ARTICLE SECTIONS
Jump To

This article is cited by 18 publications.

  1. Myongin Oh, Gabriel C. A. da Hora, Jessica M. J. Swanson. tICA-Metadynamics for Identifying Slow Dynamics in Membrane Permeation. Journal of Chemical Theory and Computation 2023, 19 (23) , 8886-8900. https://doi.org/10.1021/acs.jctc.3c00526
  2. Karen Palacio-Rodriguez, Hadrien Vroylandt, Lukas S. Stelzl, Fabio Pietrucci, Gerhard Hummer, Pilar Cossio. Transition Rates and Efficiency of Collective Variables from Time-Dependent Biased Simulations. The Journal of Physical Chemistry Letters 2022, 13 (32) , 7490-7496. https://doi.org/10.1021/acs.jpclett.2c01807
  3. Vojtěch Mlýnský, Michal Janeček, Petra Kührová, Thorben Fröhlking, Michal Otyepka, Giovanni Bussi, Pavel Banáš, Jiří Šponer. Toward Convergence in Folding Simulations of RNA Tetraloops: Comparison of Enhanced Sampling Techniques and Effects of Force Field Modifications. Journal of Chemical Theory and Computation 2022, 18 (4) , 2642-2656. https://doi.org/10.1021/acs.jctc.1c01222
  4. Lunna Li, Tommaso Casalini, Paolo Arosio, Matteo Salvalaglio. Modeling the Structure and Interactions of Intrinsically Disordered Peptides with Multiple Replica, Metadynamics-Based Sampling Methods and Force-Field Combinations. Journal of Chemical Theory and Computation 2022, 18 (3) , 1915-1928. https://doi.org/10.1021/acs.jctc.1c00889
  5. Jakub Rydzewski, Omar Valsson. Multiscale Reweighted Stochastic Embedding: Deep Learning of Collective Variables for Enhanced Sampling. The Journal of Physical Chemistry A 2021, 125 (28) , 6286-6302. https://doi.org/10.1021/acs.jpca.1c02869
  6. Shingo Ito, Yuji Sugita. Free-energy landscapes of transmembrane homodimers by bias-exchange adaptively biased molecular dynamics. Biophysical Chemistry 2024, 307 , 107190. https://doi.org/10.1016/j.bpc.2024.107190
  7. Jessica Gasparello, Marco Verona, Adriana Chilin, Roberto Gambari, Giovanni Marzaro. Assessing the interaction between hemoglobin and the receptor binding domain of SARS-CoV-2 spike protein through MARTINI coarse-grained molecular dynamics. International Journal of Biological Macromolecules 2023, 253 , 127088. https://doi.org/10.1016/j.ijbiomac.2023.127088
  8. Baltzar Stevensson, Mattias Edén. Improved reweighting protocols for variationally enhanced sampling simulations with multiple walkers. Physical Chemistry Chemical Physics 2023, 25 (33) , 22063-22078. https://doi.org/10.1039/D2CP04009C
  9. Ferdinand Horvath, Sascha Berlansky, Lena Maltan, Herwig Grabmayr, Marc Fahrner, Isabella Derler, Christoph Romanin, Thomas Renger, Heinrich Krobath. Swing‐out opening of stromal interaction molecule 1. Protein Science 2023, 32 (3) https://doi.org/10.1002/pro.4571
  10. Sandra J. Moore, Evelyne Deplazes, Ricardo L. Mancera. Influence of force field choice on the conformational landscape of rat and human islet amyloid polypeptide. Proteins: Structure, Function, and Bioinformatics 2023, 91 (3) , 338-353. https://doi.org/10.1002/prot.26432
  11. Thilini H. Gamage, Herwig Grabmayr, Ferdinand Horvath, Marc Fahrner, Doriana Misceo, William Edward Louch, Gjermund Gunnes, Helen Pullisaar, Janne Elin Reseland, Staale Petter Lyngstadaas, Asbjørn Holmgren, Silja S. Amundsen, Petr Rathner, Linda Cerofolini, Enrico Ravera, Heinrich Krobath, Claudio Luchinat, Thomas Renger, Norbert Müller, Christoph Romanin, Eirik Frengen. A single amino acid deletion in the ER Ca 2+ sensor STIM1 reverses the in vitro and in vivo effects of the Stormorken syndrome–causing R304W mutation. Science Signaling 2023, 16 (771) https://doi.org/10.1126/scisignal.add0509
  12. Johannes C. B. Dietschreit, Dennis J. Diestler, Andreas Hulm, Christian Ochsenfeld, Rafael Gómez-Bombarelli. From free-energy profiles to activation free energies. The Journal of Chemical Physics 2022, 157 (8) https://doi.org/10.1063/5.0102075
  13. Andreas Hulm, Johannes C. B. Dietschreit, Christian Ochsenfeld. Statistically optimal analysis of the extended-system adaptive biasing force (eABF) method. The Journal of Chemical Physics 2022, 157 (2) https://doi.org/10.1063/5.0095554
  14. Abhinav Gupta, Shivani Verma, Ramsha Javed, Suraj Sudhakar, Saurabh Srivastava, Nisanth N. Nair. Exploration of high dimensional free energy landscapes by a combination of temperature‐accelerated sliced sampling and parallel biasing. Journal of Computational Chemistry 2022, 43 (17) , 1186-1200. https://doi.org/10.1002/jcc.26882
  15. Katya Ahmad, Andrea Rizzi, Riccardo Capelli, Davide Mandelli, Wenping Lyu, Paolo Carloni. Enhanced-Sampling Simulations for the Estimation of Ligand Binding Kinetics: Current Status and Perspective. Frontiers in Molecular Biosciences 2022, 9 https://doi.org/10.3389/fmolb.2022.899805
  16. Harvinder Singh, Anupam Raja, Nishant Shekhar, Arushi Chauhan, Ajay Prakash, Pramod Avti, Bikash Medhi. Computational attributes of protein kinase-C gamma C2-domain & virtual screening for small molecules: elucidation from meta-dynamics simulations & free-energy calculations. Journal of Biomolecular Structure and Dynamics 2022, 10 , 1-12. https://doi.org/10.1080/07391102.2022.2077447
  17. Federica Lodesani, Maria Cristina Menziani, Shingo Urata, Alfonso Pedone. Biasing crystallization in fused silica: An assessment of optimal metadynamics parameters. The Journal of Chemical Physics 2022, 156 (19) , 194501. https://doi.org/10.1063/5.0089183
  18. Tomomi Kondo, Takehiko Sasaki, Sergi Ruiz‐Barragan, Jordi Ribas‐Ariño, Motoyuki Shiga. Refined metadynamics through canonical sampling using time‐invariant bias potential: A study of polyalcohol dehydration in hot acidic solutions. Journal of Computational Chemistry 2021, 42 (3) , 156-165. https://doi.org/10.1002/jcc.26443
  • Abstract

    Figure 1

    Figure 1. (a) FES obtained along one run after 106, 2 × 106, 3 × 106, 4 × 106 and 2 × 107 steps using balanced exponential reweighting, Tiwary reweighting and negative bias potential (red, blue and green points, respectively). The true potential (purple) is superimposed for comparison. (b–e) Time series of (b) the RMSDFES (in inset ), (c) the absolute value of the estimated free energy difference between the two minima of the double well potential, (d) the absolute value of the deviation of the estimated height of the free energy barrier from the true value of 1, and (e) the RMSDV. Same color scheme as in part a. The solid lines represent the average values of the quantities across the 72 independent runs. Shaded bands indicate the standard deviations. The positions of free energy minima and maxima are computed by fitting parabolas to the data points available in the intervals [−1.2,0.8], [−0.2,0.2], and [0.8,1.2] and taking the ordinate of the vertexes. (f–g)Time series of the position of the particle along one run where balanced exponential (f) and Tiwary (g) weights are reported according to a color scale.

    Figure 2

    Figure 2. (a) Reference FES of the alanine dipeptide as a function of the CV φ and ψ. Contour levels are plotted every 1 kbT. The black contours indicate the region within 10 kbT from the global minimum, which is used for the RMSDFES calculations. The red triangle, square and circle indicate the position of the minima C7eq and Cax and the transition state ‡, respectively. The color shades are guides for the eye. (b–f) Difference between the reference FES and those estimated using the balanced exponential, Tiwary’s, Branduardi’s, and Bonomi’s reweighting, and the negative bias, respectively, after 1.4 ns in one of the runs. Contour levels are plotted every 0.25 kbT, and the range of color shades is 4 times smaller than that in part a.

    Figure 3

    Figure 3. (a–c) Time series of the RMSDFES and the deviation from the reference ΔFC7eqCax and ΔFCax–‡, respectively, for the balanced exponential (red), Tiwary’s (blue), Branduardi’s (purple), and Bonomi’s (orange) reweighting, and negative bias (green) estimate of the FES. The data are averaged over 96 runs. Shaded bands indicate standard deviations. The black line at 1.4 ns indicates the time point where the FESs in Figure 2 have been extracted. (a, inset) shows to approximately reach a plateau for balanced exponential, Tiwary’s, and Bonomi’s reweighting schemes. (d) RMSDV between the different reweighting schemes and negative bias (same color scheme as above). (e, f) Time series of ψ when φ = −1.88 rad along one selected run. Balanced exponential (e) and Tiwary (f) weights are reported according to the color scale. ⟨w⟩ is the average of the weights along the run.

    Figure 4

    Figure 4. Kish’s effective sample size of the Tiwary’s (purple) and balanced exponential (green) reweighting method as a function of the total sample size along one standard metadynamics trajectory of the alanine dipeptide (data saved every 0.02 ps). The blue line represents the size of the ideally uncorrelated sample.

    Figure 5

    Figure 5. Block analysis of the average error on the determination of the free energy landscape of the alanine dipeptide using standard metadynamics simulations and the Tiwary’s (green) or the balanced exponential (purple) reweighting scheme. The computation is performed on a randomly selected 2 ns-long trajectory. The error is averaged over the regions of the FES where the free energy is within 10kbT of the global minimum.

    Figure 6

    Figure 6. Same as Figure 2 but for the standard metadynamics runs using only φ as biasing collective variable.

    Figure 7

    Figure 7. Same as Figure 3 but for the standard metadynamics runs using only φ as biasing collective variable.

    Figure 8

    Figure 8. Same as Figure 2 but for the well-tempered metadynamics runs.

    Figure 9

    Figure 9. Same as Figure 3 but for the well-tempered metadynamics runs. In the inset of panel a, the reaches a plateau, as expected, although this is more evident when using a larger free energy cutoff in the calculation of the RMSDFES as shown in Figure S4 in the Supporting Information.

  • References

    ARTICLE SECTIONS
    Jump To

    This article references 30 other publications.

    1. 1
      Bernardi, R. C.; Melo, M. C.; Schulten, K. Enhanced sampling techniques in molecular dynamics simulations of biological systems. Biochim. Biophys. Acta, Gen. Subj. 2015, 1850, 872877,  DOI: 10.1016/j.bbagen.2014.10.019
    2. 2
      Dror, R. O.; Dirks, R. M.; Grossman, J.; Xu, H.; Shaw, D. E. Biomolecular Simulation: A Computational Microscope for Molecular Biology. Annu. Rev. Biophys. 2012, 41, 429452,  DOI: 10.1146/annurev-biophys-042910-155245
    3. 3
      Abrams, C.; Bussi, G. Enhanced sampling in molecular dynamics using metadynamics, replica-exchange, and temperature-acceleration. Entropy 2014, 16, 163199,  DOI: 10.3390/e16010163
    4. 4
      Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23, 187199,  DOI: 10.1016/0021-9991(77)90121-8
    5. 5
      Hansmann, U. H. Parallel tempering algorithm for conformational studies of biological molecules. Chem. Phys. Lett. 1997, 281, 140150,  DOI: 10.1016/S0009-2614(97)01198-6
    6. 6
      Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 1999, 314, 141151,  DOI: 10.1016/S0009-2614(99)01123-9
    7. 7
      Sugita, Y.; Okamoto, Y. Replica-exchange multicanonical algorithm and multicanonical replica-exchange method for simulating systems with rough energy landscape. Chem. Phys. Lett. 2000, 329, 261270,  DOI: 10.1016/S0009-2614(00)00999-4
    8. 8
      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, 291318,  DOI: 10.1146/annurev.physchem.53.082301.113146
    9. 9
      Darve, E.; Rodríguez-Gómez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations. J. Chem. Phys. 2008, 128, 144120,  DOI: 10.1063/1.2829861
    10. 10
      Maragliano, L.; Vanden-Eijnden, E. A temperature accelerated method for sampling free energy and determining reaction pathways in rare events simulations. Chem. Phys. Lett. 2006, 426, 168175,  DOI: 10.1016/j.cplett.2006.05.062
    11. 11
      Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 1256212566,  DOI: 10.1073/pnas.202427399
    12. 12
      Bussi, G.; Laio, A.; Parrinello, M. Equilibrium free energies from nonequilibrium metadynamics. Phys. Rev. Lett. 2006, 96, 090601  DOI: 10.1103/PhysRevLett.96.090601
    13. 13
      Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100, 020603,  DOI: 10.1103/PhysRevLett.100.020603
    14. 14
      Bonomi, M.; Barducci, A.; Parrinello, M. Reconstructing the equilibrium Boltzmann distribution from well-tempered metadynamics. J. Comput. Chem. 2009, 30, 16151621,  DOI: 10.1002/jcc.21305
    15. 15
      Branduardi, D.; Bussi, G.; Parrinello, M. Metadynamics with adaptive Gaussians. J. Chem. Theory Comput. 2012, 8, 2247,  DOI: 10.1021/ct3002464
    16. 16
      Tiwary, P.; Parrinello, M. A Time-Independent Free Energy Estimator for Metadynamics. J. Phys. Chem. B 2015, 119, 736742,  DOI: 10.1021/jp504920s
    17. 17
      Tiana, G. Estimation of microscopic averages from metadynamics. Eur. Phys. J. B 2008, 63, 235,  DOI: 10.1140/epjb/e2008-00232-8
    18. 18
      Smiatek, J.; Heuer, A. Calculation of free energy landscapes: A histogram reweighted metadynamics approach. J. Comput. Chem. 2011, 32, 20842096,  DOI: 10.1002/jcc.21790
    19. 19
      Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185, 604613,  DOI: 10.1016/j.cpc.2013.09.018
    20. 20
      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 dynamics. Comput. Phys. Commun. 2009, 180, 1961,  DOI: 10.1016/j.cpc.2009.05.011
    21. 21
      Gimondi, I.; Tribello, G. A.; Salvalaglio, M. Building maps in collective variable space. J. Chem. Phys. 2018, 149, 104104,  DOI: 10.1063/1.5027528
    22. 22
      Dama, J. F.; Parrinello, M.; Voth, G. A. Well-tempered metadynamics converges asymptotically. Phys. Rev. Lett. 2014, 112, 240602,  DOI: 10.1103/PhysRevLett.112.240602
    23. 23
      Laio, A.; Gervasio, F. L. Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys. 2008, 71, 126601,  DOI: 10.1088/0034-4885/71/12/126601
    24. 24
      Marinelli, F.; Pietrucci, F.; Laio, A.; Piana, S. A Kinetic Model of Trp-Cage Folding from Multiple Biased Molecular Dynamics Simulations. PLoS Comput. Biol. 2009, 5, e1000452  DOI: 10.1371/journal.pcbi.1000452
    25. 25
      Laio, A.; Rodriguez-Fortea, A.; Gervasio, F. L.; Ceccarelli, M.; Parrinello, M. Assessing the Accuracy of Metadynamics †. J. Phys. Chem. B 2005, 109, 67146721,  DOI: 10.1021/jp045424k
    26. 26
      Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A message-passing parallel molecular dynamics implementation. Comput. Phys. Commun. 1995, 91, 4356,  DOI: 10.1016/0010-4655(95)00042-E
    27. 27
      Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber biomolecular simulation programs. J. Comput. Chem. 2005, 26, 16681688,  DOI: 10.1002/jcc.20290
    28. 28
      Kish, L. Survey Sampling; John Wiley & Sons, Ltd.; New York, 1965.
    29. 29
      Frenkel, D.; Smit, B. Understanding Molecular Simulations, 2nd ed.; Computational Science; Academic Press: 2002; Vol. 1.
    30. 30
      Bussi, G. Personal communication.
  • Supporting Information

    Supporting Information

    ARTICLE SECTIONS
    Jump To

    The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.9b00867.

    • Figures S1–S4 (PDF)


    Terms & Conditions

    Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system: http://pubs.acs.org/page/copyright/permissions.html.

Pair your accounts.

Export articles to Mendeley

Get article recommendations from ACS based on references in your Mendeley library.

Pair your accounts.

Export articles to Mendeley

Get article recommendations from ACS based on references in your Mendeley library.

You’ve supercharged your research process with ACS and Mendeley!

STEP 1:
Click to create an ACS ID

Please note: If you switch to a different device, you may be asked to login again with only your ACS ID.

Please note: If you switch to a different device, you may be asked to login again with only your ACS ID.

Please note: If you switch to a different device, you may be asked to login again with only your ACS ID.

MENDELEY PAIRING EXPIRED
Your Mendeley pairing has expired. Please reconnect