###### Abstract

We review a selection of methods for performing enhanced sampling in molecular
dynamics simulations. We consider methods based on collective variable biasing
and on tempering, and offer both historical and contemporary perspectives. In
collective-variable biasing, we first discuss methods stemming from
thermodynamic integration that use mean force biasing, including the adaptive
biasing force algorithm and temperature acceleration. We then turn to methods
that use bias potentials, including umbrella sampling and metadynamics. We
next consider parallel tempering and replica-exchange methods. We conclude
with a brief presentation of some combination methods.

###### keywords:

collective variables, free energy, blue-moon sampling, adaptive-biasing force
algorithm, temperature-acceleration, umbrella sampling, metadynamics

10.3390/—— xx Received: xx / Accepted: xx / Published: xx Enhanced sampling in
molecular dynamics using metadynamics, replica-exchange, and temperature-
acceleration Cameron Abrams 1,⋆ and Giovanni Bussi 2<EMAIL_ADDRESS>tel
+1-215-895-2231

##  1 Introduction

The purpose of molecular dynamics (MD) is to compute the positions and
velocities of a set of interacting atoms at the present time instant given
these quantities one time increment in the past. Uniform sampling from the
discrete trajectories one can generate using MD has long been seen as
synonymous with sampling from a statistical-mechanical ensemble; this just
expresses our collective wish that the ergodic hypothesis holds at finite
times. Unfortunately, most MD trajectories are not ergodic and leave many
relevant regions of configuration space unexplored. This stems from the
separation of high-probability “metastable” regions by low-probability
“transition” regions and the inherent difficulty of sampling a
3$N$-dimensional space by embedding into it a one-dimensional dynamical
trajectory.

This review concerns a selection of methods to use MD simulation to enhance
the sampling of configuration space. A central concern with any enhanced
sampling method is guaranteeing that the statistical weights of the samples
generated are known and correct (or at least correctable) while simultaneously
ensuring that as much of the relevant regions of configuration space are
sampled. Because of the tight relationship between probability and free
energy, many of these methods are known as “free-energy” methods. To be sure,
there are a large number of excellent reviews of free-energy methods in the
literature (e.g., Kollman (1993); Trzesniak et al. (2007); Vanden-Eijnden
(2009); Dellago and Bolhuis (2009); Christ et al. (2010)). The present review
is in no way intended to be as comprehensive as these; as the title indicates,
we will mostly focus on enhanced sampling methods of three flavors: tempering,
metadynamics, and temperature-acceleration. Along the way, we will point out
important related methods, but in the interest of brevity we will not spend
much time explaining these. The methods we have chosen to focus on reflect our
own preferences to some extent, but they also represent popular and growing
classes of methods that find ever more use in biomolecular simulations and
beyond.

We divide our review into three main sections. In the first, we discuss
enhanced sampling approaches that rely on collective variable biasing. These
include the historically important methods of thermodynamic integration and
umbrella sampling, and we pay particular attention to the more recent
approaches of the adaptive-biasing force algorithm, temperature-acceleration,
and metadynamics. In the second section, we discuss approaches based on
tempering, which is dominated by a discussion of the parallel
tempering/replica exchange approaches. In the third section, we briefly
present some relatively new methods derived from either collective-variable-
based or tempering-based approaches, or their combinations.

##  2 Approaches Based on Collective-Variable Biasing

###  2.1 Background: Collective Variables and Free Energy

For our purposes, the term “collective variable” or CV refers to any
multidimensional function ${\boldsymbol{\theta}}$ of 3$N$-dimensional atomic
configuration ${\boldsymbol{x}}\equiv\left(x_{i}|i=1\dots 3N\right)$. The
functions $\theta_{1}({\boldsymbol{x}})$,
$\theta_{2}({\boldsymbol{x}})$,$\dots$,$\theta_{M}({\boldsymbol{x}})$ map
configuration ${\boldsymbol{x}}$ onto an $M$-dimensional CV space
${\boldsymbol{z}}\equiv\left(z_{j}|j=1\dots M\right)$, where usually $M\ll
3N$. At equilibrium, the probability of observing the system at CV-point
${\boldsymbol{z}}$ is the weight of all configurations ${\boldsymbol{x}}$
which map to ${\boldsymbol{z}}$:

$P({\boldsymbol{z}})=\left<\delta[{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}]\right>,$
(1)

The Dirac delta function picks out only those configurations for which the CV
${\boldsymbol{\theta}}({\boldsymbol{x}})$ is ${\boldsymbol{z}}$, and
$\left<\cdot\right>$ denotes averaging its argument over the equilibrium
probability distribution of ${\boldsymbol{x}}$. The probability can be
expressed as a free energy:

$F({\boldsymbol{z}})=-k_{B}T\ln\left<\delta[{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}]\right>.$
(2)

Here, $k_{B}$ is Boltzmann’s constant and $T$ is temperature.

Local minima in $F$ are metastable equilibrium states. $F$ also measures the
energetic cost of a maximally efficient (i.e., reversible) transition from one
region of CV space to another. If, for example, we choose a CV space such that
two well-separated regions define two important allosteric states of a given
protein, we could perform a free-energy calculation to estimate the change in
free energy required to realize the conformational transition. Indeed, the
promise of being able to observe with atomic detail the transition states
along some pathway connecting two distinct states of a biomacromolecule is
strong motivation for exploring these transitions with CV’s.

Given the limitations of standard MD, how does one “discover” such states in a
proposed CV space? A perfectly ergodic (infinitely long) MD trajectory would
visit these minima much more frequently than it would the intervening spaces,
allowing one to tally how often each point in CV space is visited; normalizing
this histogram into a probability $P({\boldsymbol{z}})$ would be the most
straightforward way to compute $F$ via Eq. 2. In all too many actual cases, MD
trajectories remain close to only one minimum (the one closest to the initial
state of the simulation) and only very rarely, if ever, visit others. In the
CV sense, we therefore speak of standard MD simulations failing to overcome
barriers in free energy. “Enhanced sampling” in this context refers then to
methods by which free-energy barriers in a chosen CV space are surmounted to
allow as broad as possible an extent of CV space to be explored and
statistically characterized with limited computational resources.

In this section, we focus on methods of enhanced sampling of CV’s based on MD
simulations that are directly biased on those CV’s; that is, we focus on
methods in which an investigator must identify the CV’s of interest as an
input to the calculation. We have chosen to limit discussion to two broad
classes of biasing: those whose objective is direct computation of the
gradient of the free energy $(\partial F/\partial{\boldsymbol{z}})$ at local
points throughout CV space, and those in which non-Boltzmann sampling with
bias potentials is used to force exploration of otherwise hard-to-visit
regions of CV space. The canonical methods in these two classes are
thermodynamic integration and umbrella sampling, respectively, and a
discussion of these two methods sets the stage for discussion of three
relatively modern variants: the Adaptive-Biasing Force Algorithm Darve et al.
(2008), Temperature-Accelerated MD Maragliano and Vanden-Eijnden (2006) and
Metadynamics Laio and Parrinello (2002).

###  2.2 Gradient Methods: Blue-Moon Sampling, Adaptive-Biasing Force
Algorithm, and Temperature-Accelerated MD

####  2.2.1 Overview: Thermodynamic Integration

Naively, one way to have an MD system visit a hard-to-reach point
${\boldsymbol{z}}$ in CV space is simply to create a realization of the
configuration ${\boldsymbol{x}}$ at that point (i.e., such that
${\boldsymbol{\theta}}({\boldsymbol{x}})={\boldsymbol{z}}$). This is an
inverse problem, since the number of degrees of freedom in ${\boldsymbol{x}}$
is usually much larger than in ${\boldsymbol{z}}$. One way to perform this
inversion is by introducing external forces that guide the configuration to
the desired point from some easy-to-create initial state; both targeted MD
Schlitter et al. (1993) and steered MD Grubmüller et al. (1996) are ways to do
this. Of course, one would like MD to explore CV space in the vicinity of
${\boldsymbol{z}}$, so after creating the configuration ${\boldsymbol{x}}$,
one would just let it run. Unfortunately, this would likely result in the
system drifting away from ${\boldsymbol{z}}$ rather quickly, and there would
be no way from such calculations to estimate the likelihood of observing an
unbiased long MD simulation visit ${\boldsymbol{z}}$. But there is information
in the fact that the system drifts away; if one knows on average which
direction and how strongly the system would like to move if initialized at
${\boldsymbol{z}}$, this would be a measure of negative gradient of the free
energy, $-(\partial F/\partial{\boldsymbol{z}})$, or the “mean force”. We have
then a glimpse of a three-step method to compute $F$ (i.e., the statistics of
CV’s) over a meaningfully broad extent of CV space:

  1. 1.

visit a select number of local points in that space, and at each one,

  2. 2.

compute the mean force, then

  3. 3.

use numerical integration to reconstruct $F$ from these local mean forces;
formally expressed as

$F({\boldsymbol{z}})-F({\boldsymbol{z}}_{0})=\int_{{\boldsymbol{z}}_{0}}^{{\boldsymbol{z}}}\left(\frac{\partial
F}{\partial{\boldsymbol{z}}}\right)d{\boldsymbol{z}}$ (3)

Inspired by Kirkwood’s original suggestion involving switching parameters
Kirkwood (1935), such an approach is generally referred to as “thermodynamic
integration” or TI. TI allows us to reconstruct the statistical weights of any
point in CV space by accumulating information on the gradients of free energy
at selected points.

####  2.2.2 Blue-Moon Sampling

The discussion so far leaves open the correct way to compute the local free-
energy gradients. A gradient is a local quantity, so a natural choice is to
compute it from an MD simulation localized at a point in CV space by a
constraint. Consider a long MD simulation with a holonomic constraint fixing
the system at the point ${\boldsymbol{z}}$. Uniform samples from this
constrained trajectory ${\boldsymbol{x}}(t)$ then represent an ensemble at
fixed ${\boldsymbol{z}}$ over which the averaging needed to convert gradients
in potential energy to gradients in free energy could be done. However, this
constrained ensemble has the undesired property that the velocities
$\dot{\boldsymbol{\theta}}({\boldsymbol{x}})$ are zero. This is a bit
problematic because virtually none of the samples plucked from a long
unconstrained MD simulation (as is implied by Eq. 1), would have
$\dot{\boldsymbol{\theta}}=0$, and $\dot{\boldsymbol{\theta}}=0$ acts as a set
of $M$ unphysical constraints on the system velocities $\dot{\boldsymbol{x}}$,
since $\dot{\theta}_{j}=\sum_{i}(\partial\theta_{j}/\partial
x_{i})\dot{x}_{i}$. Probably the best-known example of a method to correct for
this bias is the so-called “blue-moon” sampling method Carter et al. (1989);
Sprik and Ciccotti (1998); Ciccotti and Ferrario (2004); Ciccotti et al.
(2005) or the constrained ensemble method den Otter and Briels (1998);
Schlitter and Klähn (2003). The essence of the method is a decomposition of
free energy gradients into components along the CV gradients and thermal
components orthogonal to them:

$\frac{\partial F}{\partial
z_{j}}=\left<{\boldsymbol{b}}_{j}({\boldsymbol{x}})\cdot\nabla
V({\boldsymbol{x}})-k_{B}T\nabla\cdot{\boldsymbol{b}}_{j}({\boldsymbol{x}})\right>_{{\boldsymbol{\theta}}({\boldsymbol{x}})={\boldsymbol{z}}}$
(4)

where
$\left<\cdot\right>_{{\boldsymbol{\theta}}({\boldsymbol{x}})={\boldsymbol{z}}}$
denotes averaging across samples drawn uniformly from the MD simulation
constrained at ${{\boldsymbol{\theta}}({\boldsymbol{x}})={\boldsymbol{z}}}$,
and the ${\boldsymbol{b}}_{j}({\boldsymbol{x}})$ is the vector field
orthogonal to the gradients of every component $k$ of ${\boldsymbol{\theta}}$
for $k\neq j$:

${\boldsymbol{b}}_{j}({\boldsymbol{x}})\cdot\nabla\theta_{k}({\boldsymbol{x}})=\delta_{jk}$
(5)

where $\delta_{jk}$ is the Kroenecker delta. (For brevity, we have omitted the
consideration of holonomic constraints other than that on the CV; the reader
is referred to the paper by Ciccotti et al. for details Ciccotti et al.
(2005).) The vector fields ${\boldsymbol{b}}_{j}$ for each $\theta_{j}$ can be
constructed by orthogonalization. The first term in the angle brackets in Eq.
4 implements the chain rule one needs to account for how energy $V$ changes
with ${\boldsymbol{z}}$ through all the ways ${\boldsymbol{z}}$ can change
with ${\boldsymbol{x}}$. The second term corrects for the thermal bias imposed
by the constraint.

Although nowhere near exhaustive, below is a listing of common types of
problems to which blue-moon sampling has been applied with some representative
examples:

  1. 1.

sampling conformations of small flexible molecules and peptides Depaepe et al.
(1993); Zhao et al. (2008); Kim et al. (2009)

  2. 2.

environmental effects on covalent bond formation/breaking (usually in
combination with ab initio MD) Hytha et al. (2001); Fois et al. (2004); Ivanov
and Klein (2005); Stubbs and Marx (2005); Trinh et al. (2009); Liu et al.
(2010); Bucko and Hafner (2010)

  3. 3.

solvation and non-covalent binding of small molecules in solvent Paci and
Marchi (1994); Sa et al. (2006); Mugnai et al. (2007); Chunsrivirot and Trout
(2011); Sato et al. (2012)

  4. 4.

protein dimerization Sergi et al. (2002); Maragliano et al. (2004)

####  2.2.3 The Adaptive Biasing Force Algorithm

The blue-moon approach requires multiple independent constrained MD
simulations to cover the region of CV space in which one wants internal
statistics. The care taken in choosing these quadrature points can often
dictate the accuracy of the resulting free energy reconstruction. It is
therefore sometimes advantageous to consider ways to avoid having to choose
such points ahead of time, and adaptive methods attempt to address this
problem. One example is the adaptive-biasing force (ABF) algorithm of Darve et
al. Darve and Pohorille (2001); Darve et al. (2008) The essence of ABF is two-
fold: (1) recognition that external bias forces of the form
$\nabla_{\boldsymbol{x}}\theta_{j}\left(\partial F/\partial z_{j}\right)$ for
$j=1\dots M$ exactly oppose mean forces and should lead to more uniform
sampling of CV space, and (2) that these bias forces can be converged upon
adaptively during a single unconstrained MD simulation.

The first of those two ideas is motivated by the fact that “forces” that keep
normal MD simulations effectively confined to free energy minima are mean
forces on the collective variables projected onto the atomic coordinates, and
balancing those forces against their exact opposite should allow for thermal
motion to take the system out of those minima. The second idea is a bit more
subtle; after all, in a running MD simulation with no CV constraints, the
constrained ensemble expression for the mean force (Eq. 4) does not directly
apply, because a constrained ensemble is not what is being sampled. However,
Darve et al. showed how to relate these ensembles so that the samples
generated in the MD simulation could be used to build mean forces Darve and
Pohorille (2001). Further, they showed using a clever choice of the fields of
Eq. 4 an equivalence between ($i$) the spatial gradients needed to computed
forces, and ($ii$) time-derivatives of the CV’s Darve et al. (2008):

$\frac{\partial F}{\partial
z_{i}}=-k_{B}T\left<\frac{d}{dt}\left(M_{\theta}\frac{d\theta_{i}}{dt}\right)\right>_{{\boldsymbol{\theta}}={\boldsymbol{z}}}$
(6)

where $M_{\theta}$ is the transformed mass matrix given by

$M_{\theta}^{-1}=J_{\theta}M^{-1}J_{\theta}$ (7)

where $J_{\theta}$ is the $M\times 3N$ matrix with elements
$\partial\theta_{i}/\partial x_{j}$ ($i=1\dots M$, $j=1\dots 3N$), and $M$ is
the diagonal matrix of atomic masses. Eq. 7 is the result of a particular
choice for the fields ${\boldsymbol{b}}_{j}({\boldsymbol{x}})$. This
reformulation of the instantaneous mean forces computed on-the-fly makes ABF
exceptionally easy to implement in most modern MD packages. Darve et al.
present a clear demonstration of the ABF algorithm in a  pseudocode Darve et
al. (2008) that attests to this fact.

ABF has found rather wide application in CV-based free energy calculations in
recent years. Below is a  representative sample of some types of problems
subjected to ABF calculations in the recent literature:

  1. 1.

Peptide backbone angle sampling Fogolari et al. (2011); Faller et al. (2013);

  2. 2.

Nucleoside Wei and Pohorile (2011), protein Vivcharuk and Kaznessis (2011) and
fullerene Kraszewski et al. (2011, 2012) insertion into a lipid bilayer;

  3. 3.

Interactions of small molecules with polymers in water Liu et al. (2010);
Caballero et al. (2013);

  4. 4.

Molecule/ion transport through protein complexes Wilson et al. (2011); Cheng
and Coalson (2012); Wang et al. (2012); Tillman et al. (2013) and DNA
superstructures Akhshi et al. (2012);

  5. 5.

Calculation of octanol-water partition coefficients Kamath et al. (2012);
Bhatnagar et al. (2012);

  6. 6.

Large-scale protein conformational changes Wereszczynski and McCammon (2012);

  7. 7.

Protein-nanotube Jana and Sengupta (2012) and nanotube-nanotube Uddin et al.
(2010) association.

####  2.2.4 Temperature-Accelerated Molecular Dynamics

Both blue-moon sampling and ABF are based on statistics in the constrained
ensemble. However, estimation of mean forces need not only use this ensemble.
One can instead relax the constraint and work with a “mollified” version of
the free energy:

$F_{\kappa}({\boldsymbol{z}})=-k_{B}T\ln\left<\delta_{\kappa}\left[{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}\right]\right>$
(8)

where $\delta_{\kappa}$ refers to the Gaussian (or “mollified delta
function”):

$\delta_{\kappa}=\sqrt{\frac{\beta\kappa}{2\pi}}\exp\left[-\frac{1}{2}\beta\kappa\left|{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}\right|^{2}\right],$
(9)

where $\beta$ is just shorthand for $1/k_{B}T$. Since
$\lim_{\beta\kappa\rightarrow\infty}\delta_{\kappa}=\delta$, we know that
$\lim_{\beta\kappa\rightarrow\infty}F_{\kappa}=F$. One way to view this
Gaussian is that it “smoothes out” the true free energy to a tunable degree;
the factor $1/\sqrt{\beta\kappa}$ is a length-scale in CV space below which
details are smeared.

Because the Gaussian has continuous gradients, it can be used directly in an
MD simulation. Suppose we have a CV space
${\boldsymbol{\theta}}({\boldsymbol{x}})$, and we extend our MD system to
include variables ${\boldsymbol{z}}$ such that the combined set
$({\boldsymbol{x}},{\boldsymbol{z}})$ obeys the following extended potential:

$U({\boldsymbol{x}},{\boldsymbol{z}})=V({\boldsymbol{x}})+\sum_{j=1}^{M}\frac{1}{2}\kappa\left|\theta_{j}({\boldsymbol{x}})-z_{j}\right|^{2}$
(10)

where $V({\boldsymbol{x}})$ is the interatomic potential, and $\kappa$ is a
constant. Clearly, if we fix ${\boldsymbol{z}}$, then the resulting free
energy is to within an additive constant the mollified free energy of Eq. 8.
(The additive constant is related to the prefactor of the mollified delta
function and has nothing to do with the number of CV’s.) Further, we can
directly express the gradient of this mollified free energy with respect to
${\boldsymbol{z}}$: Kaestner (2009)

$\nabla_{\boldsymbol{z}}F_{\kappa}=-\left<\kappa\left[{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}\right]\right>$
(11)

This suggests that, instead of using constrained ensemble MD to accumulate
mean forces, we could work in the restrained ensemble and get very good
approximations to the mean force. By “restrained”, we refer to the fact that
the term giving rise to the mollified delta function in the configurational
integral is essentially a harmonic restraining potential with a “spring
constant” $\kappa$. In this restrained-ensemble approach, no velocities are
held fixed, and the larger we choose $\kappa$ the more closely we can
approximate the true free energy. Notice however that large values of $\kappa$
could lead to numerical instabilities in integrating equations of motion, and
a balance should be found. (In practice, we have found that for CV’s with
dimensions of length, values of $\kappa$ less than about 1,000 kcal/mol/Å2 can
be stably handled, and values of around 100 kcal/mol/Å2 are typically
adequate.)

Temperature-accelerated MD (TAMD) Maragliano and Vanden-Eijnden (2006) takes
advantage of the restrained-ensemble approach to directly evolve the variables
${\boldsymbol{z}}$ in such a way to accelerate the sampling of CV space.
First, consider how the atomic variables ${\boldsymbol{x}}$ evolve under the
extended potential (assuming Langevin dynamics):

$m_{i}\ddot{x}_{i}=-\frac{\partial V({\boldsymbol{x}})}{\partial
x_{i}}-\kappa\sum_{j=1}^{m}\left[\theta_{j}({\boldsymbol{x}})-z_{j}\right]\frac{\partial\theta_{j}({\boldsymbol{x}})}{\partial
x_{i}}-\gamma m_{i}\dot{x_{i}}+\eta_{i}(t;\beta)$ (12)

Here, $m_{i}$ is the mass of $x_{i}$, $\gamma$ is the friction coefficient for
the Langevin thermostat, and ${\boldsymbol{\eta}}$ is the thermostat white
noise satisfying the fluctuation-dissipation theorem at physical temperature
$\beta^{-1}$:

$\left<\eta_{i}(t;\beta)\eta_{j}(t^{\prime};\beta)\right>=\beta^{-1}\gamma
m_{i}\delta_{ij}\delta(t-t^{\prime})$ (13)

Key to TAMD is that the ${\boldsymbol{z}}$ are treated as slow variables that
evolve according to their own equations of motion, which here we take as
diffusive (though other choices are possible Maragliano and Vanden-Eijnden
(2006)):

$\bar{\gamma}\bar{m}_{j}\dot{z}_{j}=\kappa\left[\theta_{j}({\boldsymbol{x}})-z_{j}\right]+\xi_{j}(t;\bar{\beta}).$
(14)

Here, $\bar{\gamma}$ is a fictitious friction, $\bar{m}_{j}$ is a mass, and
the first term on the right-hand side represents the instantaneous force on
variable $z_{j}$, and the second term represents thermal noise at the
fictitious thermal energy $\bar{\beta}^{-1}\not=\beta^{-1}$.

The advantage of TAMD is that if (1) $\bar{\gamma}$ is chosen sufficiently
large so as to guarantee that the slow variables indeed evolve slowly relative
to the fundamental variables, and (2) $\kappa$ is sufficiently large such that
${\boldsymbol{\theta}}({\boldsymbol{x}}(t))\approx{\boldsymbol{z}}(t)$ at any
given time, then the force acting on ${\boldsymbol{z}}$ is approximately equal
to minus the gradient of the free energy (Eq. 11) Maragliano and Vanden-
Eijnden (2006). This is because the MD integration repeatedly samples
$\kappa\left[{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}\right]$
for an essentially fixed (but actually very slowly moving) ${\boldsymbol{z}}$,
so ${\boldsymbol{z}}$ evolution effectively feels these samples as a mean
force. In other words, the dynamics of ${\boldsymbol{z}}(t)$ is effectively

$\bar{\gamma}\bar{m}_{j}\dot{z}_{j}=-\frac{\partial
F({\boldsymbol{z}})}{\partial z_{j}}+\xi_{j}(t;\bar{\beta}).$ (15)

This shows that the ${\boldsymbol{z}}$-dynamics describes an equilibrium
constant-temperature ensemble at fictitious temperature $\bar{\beta}^{-1}$
acted on by the “potential” $F({\boldsymbol{z}})$, which is the free energy
evaluated at the physical temperature $\beta^{-1}$. That is, under TAMD,
${\boldsymbol{z}}$ conforms to a probability distribution of the form
$\exp\left[-\bar{\beta}F({\boldsymbol{z}};\beta)\right]$, whereas under normal
MD it would conform to $\exp\left[-\beta F({\boldsymbol{z}};\beta)\right]$.
The all-atom MD simulation (at $\beta$) simply serves to approximate the local
gradients of $F({\boldsymbol{z}})$. Sampling is enhanced by taking
$\bar{\beta}^{-1}>\beta^{-1}$, which has the effect of attenuating the
ruggedness of $F$. TAMD therefore can accelerate a trajectory
${\boldsymbol{z}}(t)$ through CV space by increasing the likelihood of
visiting points with relatively low physical Boltzmann factors. This borrows
directly from the main idea of adiabatic free-energy dynamics Rosso and
Tuckerman (2002) (AFED), in that one deliberately makes some variables hot (to
overcome barriers) but slow (to keep them adiabatically separated from all
other variables). In TAMD, however, the use of the mollified free energy means
no cumbersome variable transformations are required. (The authors of AFED
refer to TAMD as “driven”-AFED, or d-AFED Abrams and Tuckerman (2008).) It is
also worth mentioning in this review that TAMD borrows heavily from an early
version of metadynamics Iannuzzi et al. (2003), which was formulated as a way
to evolve the auxiliary variables ${\boldsymbol{z}}$ on a mollified free
energy. However, unlike metadynamics (which we discuss below in Sec. 2.3.3),
there is no history-dependent bias in TAMD.

Unlike TI, ABF, and the methods of umbrella sampling and metadynamics
discussed in the next section, TAMD is not a method for direct calculation of
the free energy. Rather, it is a way to overcome free energy barriers in a
chosen CV space quickly without visiting irrelevant regions of CV space.
(However, we discuss briefly a method in Sec. 4.2.2 in which TAMD gradients
are used in a spirit similar to ABF to reconstruct a free energy.) That is, we
consider TAMD a way to efficiently explore relevant regions CV space that are
practically inaccessible to standard MD simulation. It is also worth pointing
out that, unlike ABF, TAMD does not operate by opposing the natural gradients
in free energy, but rather by using them to guide accelerated sampling. ABF
can only use forces in locations in CV space the trajectory has visited, which
means nothing opposes the trajectory going to regions of very high free
energy. However, under TAMD, an acceleration of $\bar{\beta}^{-1}$= 6 kcal/mol
on the CV’s will greatly accelerate transitions over barriers of 6-12
kcal/mol, but will still not (in theory) accelerate excursions to regions
requiring climbs of hundreds of kcal/mol. TAMD and ABF have in common the
ability to handle rather high-dimensional CV’s.

Although it was presented theoretically in 2006 Maragliano and Vanden-Eijnden
(2006), TAMD was not applied directly to large-scale MD until much later
Abrams and Vanden-Eijnden (2010). Since then, there has been growing interest
in using TAMD in a variety of applications requiring enhanced sampling:

  1. 1.

TAMD-enhanced flexible fitting of all-atom protein and RNA models into low-
resolution electron microscopy density maps Vashisth et al. (2012, 2013);

  2. 2.

Large-scale (interdomain) protein conformational sampling Abrams and Vanden-
Eijnden (2010); Vashisth and Brooks (2012); Hu et al. (2012);

  3. 3.

Loop conformational sampling in proteins Vashisth and Abrams (2012);

  4. 4.

Mapping of diffusion pathways for small molecules in globular proteins
Maragliano et al. (2010); Lapelosa and Abrams (2013);

  5. 5.

Vacancy diffusion Geslin et al. (2013);

  6. 6.

Conformational sampling and packing in dense polymer systems Lucid et al.
(2013).

Finally, we mention briefly that TAMD can be used as a quick way to generate
trajectories from which samples can be drawn for subsequent mean-force
estimation for later reconstruction of a multidimensional free energy; this is
the essence of the single-sweep method Maragliano and Vanden-Eijnden (2008),
which is an efficient means of computing multidimensional free energies.
Rather than using straight numerical TI, single sweep posits the free energy
as a basis function expansion and uses standard optimization methods to find
the expansion coefficients that best reproduce the measured mean forces.
Single-sweep has been used to map diffusion pathways of CO and H2O in
myoglobin Maragliano et al. (2010); Lapelosa and Abrams (2013).

###  2.3 Bias Potential Methods: Umbrella Sampling and Metadynamics

####  2.3.1 Overview: Non-Boltzmann Sampling

In the previous section, we considered methods that achieve enhanced sampling
by using mean forces: in TI, these are integrated to reconstruct a free
energy; in ABF, these are built on-the-fly to drive uniform CV sampling; and
in TAMD, these are used on-the-fly to guide accelerated evolution of CV’s. In
this section, we consider methods that achieve enhanced sampling by means of
controlled bias potentials. As a class, we refer to these as non-Boltzmann
sampling methods.

Non-Boltzmann sampling is generally a way to derive statistics on a system
whose energetics differ from the energetics used to perform the sampling.
Imagine we have an MD system with bare interatomic potential
$V({\boldsymbol{x}})$, and we add a bias $\Delta V({\boldsymbol{x}})$ to
arrive at a biased total potential:

$V_{b}({\boldsymbol{x}})=V({\boldsymbol{x}})+\Delta V({\boldsymbol{x}})$ (16)

The statistics on the CV’s on this biased potential are then given as

$\displaystyle P_{b}({\boldsymbol{z}})$
$\displaystyle=\frac{\displaystyle\int\\!d{\boldsymbol{x}}\ e^{-\beta
V_{0}({\boldsymbol{x}})}e^{-\beta\Delta
V({\boldsymbol{x}})}\delta\left[{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}\right]}{\displaystyle\int\\!d{\boldsymbol{x}}\
e^{-\beta V_{0}({\boldsymbol{x}})}e^{-\beta\Delta V({\boldsymbol{x}})}}$
$\displaystyle=\frac{\displaystyle\displaystyle\int\\!d{\boldsymbol{x}}\
e^{-\beta V_{0}({\boldsymbol{x}})}e^{-\beta\Delta
V}\delta\left[{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}\right]}{\displaystyle\int\\!d{\boldsymbol{x}}\
e^{-\beta V_{0}({\boldsymbol{x}})}}\frac{\int\\!d{\boldsymbol{x}}\ e^{-\beta
V_{0}({\boldsymbol{x}})}}{\displaystyle\int\\!d{\boldsymbol{x}}\ e^{-\beta
V_{0}({\boldsymbol{x}})}e^{-\beta\Delta V({\boldsymbol{x}})}}$
$\displaystyle=\frac{\left<e^{-\beta\Delta
V({\boldsymbol{x}})}\delta\left[{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}\right]\right>}{\left<e^{-\beta\Delta
V({\boldsymbol{x}})}\right>}$ (17)

where $\left<\cdot\right>$ denotes ensemble averaging on the unbiased
potential $V({\boldsymbol{x}})$. Further, if we take the bias potential
$\Delta V$ to be explicitly a function only of the CV’s
${\boldsymbol{\theta}}$, then it becomes invariant in the averaging of the
numerator thanks to the delta function, and we have

$P_{b}({\boldsymbol{x}})=\frac{e^{-\beta\Delta
V({\boldsymbol{z}})}\left<\delta\left[{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}\right]\right>}{\left<e^{-\beta\Delta
V\left[\theta({\boldsymbol{x}})\right]}\right>}$ (18)

Finally, since the unbiased statistics are
$P({\boldsymbol{z}})=\left<\delta\left[{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}\right]\right>$,
we arrive at

$P({\boldsymbol{z}})=P_{b}({\boldsymbol{z}})e^{\beta\Delta
V({\boldsymbol{z}})}\left<e^{-\beta\Delta
V\left[\theta({\boldsymbol{x}})\right]}\right>$ (19)

Taking  samples from an ergodic MD simulation on  the biased potential
$V_{b}$, Eq. 19 provides the recipe for reconstructing the statistics the CV’s
would present were they generated using the unbiased potential $V$. However,
the probability $P({\boldsymbol{z}})$ is  implicit in this equation, because

$\left<e^{-\beta\Delta V}\right>=\int
d{\boldsymbol{z}}P({\boldsymbol{z}})e^{-\beta\Delta
V\left[{\boldsymbol{\theta}}({\boldsymbol{x}})\right]}$ (20)

This is not really a problem, since we can treat $\left<e^{-\beta\Delta
V}\right>$ as a constant we can get from normalizing
$P_{b}({\boldsymbol{z}})e^{\beta\Delta V({\boldsymbol{z}})}$.

How does one choose $\Delta V$ so as to enhance the sampling of CV space?
Evidently, from the standpoint of non-Boltzmann sampling, the closer the bias
potential is to the negative free energy $-F({\boldsymbol{z}})$, the more
uniform the sampling of CV space will be. To wit: if $\Delta
V\left[{\boldsymbol{\theta}}({\boldsymbol{x}})\right]=-F\left[{\boldsymbol{\theta}}({\boldsymbol{x}})\right]$,
then $e^{\beta\Delta V({\boldsymbol{z}})}=e^{-\beta
F({\boldsymbol{z}})}=P({\boldsymbol{z}})$, and Eq. 19 can be inverted for
$P_{b}$ to yield

$P_{b}({\boldsymbol{z}})=\frac{1}{\left<e^{\beta
F({\boldsymbol{z}})}\right>}=\frac{1}{\displaystyle\int
d{\boldsymbol{z}}P({\boldsymbol{z}})e^{\beta
F({\boldsymbol{z}})}}=\frac{1}{\displaystyle\int d{\boldsymbol{z}}e^{-\beta
F}e^{\beta F}}=\frac{1}{\displaystyle\int d{\boldsymbol{z}}}$ (21)

So we see that taking the bias potential to be the negative free energy makes
all states ${\boldsymbol{z}}$ in CV space equiprobable. This is indeed the
limit to which ABF strives by applying negative mean forces, for example Darve
et al. (2008).

We usually do not know the free energy ahead of time; if we did, we would
already know the statistics of CV space and no enhanced sampling would be
necessary. Moreover, perfectly uniform sampling of the entire CV space is
usually far from necessary, since most CV spaces have many irrelevant regions
that should be ignored. And in reference to the mean-force methods of the last
section, uniform sampling is likely not necessary to achieve accurate mean
force values; how good an estimate of $\nabla F$ is at some point
${\boldsymbol{z}}_{0}$ should not depend on how well we sampled at some other
point ${\boldsymbol{z}}_{1}$. Yet achieving uniform sampling is an
idealization since, if we do, this means we know the free energy. We now
consider two other biasing methods that aim for this ideal, either in
relatively small regions of CV space using fixed biases, or over broader
extents using adaptive biases.

####  2.3.2 Umbrella Sampling

Umbrella sampling is the standard way of using non-Boltzmann sampling to
overcome free energy barriers. In its debut Torrie and Valleau (1977),
umbrella sampling used a function $w({\boldsymbol{x}})$ that weights hard-to-
sample configurations, equivalent to adding a bias potential of the form

$\Delta V({\boldsymbol{x}})=-k_{B}T\ln w({\boldsymbol{x}}).$ (22)

$w$ is found by trial-and-error such that configurations that are easy to
sample on the unbiased potential are still easy to sample; that is, $w$ acts
like an “umbrella” covering both the easy- and hard-to-sample regions of
configuration space. Nearly always, $w$ is an explicit function of the CV’s,
$w({\boldsymbol{x}})=W[{\boldsymbol{\theta}}({\boldsymbol{x}})]$.

Coming up with the umbrella potential that would enable exploration of CV
space with a single umbrella sampling simulation that takes the system far
from its initial point is not straightforward. Akin to TI, it is therefore
advantageous to combine results from several independent trajectories, each
with its own umbrella potential that localizes it to a small volume of CV
space that overlaps with nearby volumes. The most popular way to combine the
statistics of such a set of independent umbrella sampling runs is the
weighted-histogram analysis method (WHAM) Kumar et al. (1992).

To compute statistics of CV space using WHAM, one first chooses the points in
CV space that define the little local neighborhoods, or “windows” to be
sampled and chooses the bias potential used to localize the sampling. Not
knowing how the free energy changes in CV space makes the first task somewhat
challenging, since more densely packed windows are preferred in regions where
the free energy changes rapidly; however, since the calculations are
independent, more can be added later if needed. A convenient choice for the
bias potential is a simple harmonic spring that tethers the trajectory to a
reference point ${\boldsymbol{z}}_{i}$ in CV space:

$\Delta
V_{i}({\boldsymbol{x}})=\frac{1}{2}\kappa\left|{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}_{i}\right|^{2}$
(23)

which means the dynamics of the atomic variables ${\boldsymbol{x}}$ are
identical to Eq. 12 at fixed ${\boldsymbol{z}}={\boldsymbol{z}}_{i}$. The
points $\left\\{{\boldsymbol{z}}_{i}\right\\}$ and the value of $\kappa$
(which may be point-dependent) must be chosen such that
${\boldsymbol{\theta}}\left[{\boldsymbol{x}}(t)\right]$ from any one window’s
trajectory makes excursions into the window of each of its nearest neighbors
in CV space.

Each window-restrained trajectory is directly histogrammed to yield apparent
(i.e., biased) statistics on ${\boldsymbol{\theta}}$; let us call the biased
probability in the $i$th window $P_{b,i}({\boldsymbol{z}})$. Eq. 19 again
gives the recipe to reconstruct the unbiased statistics
$P_{i}({\boldsymbol{z}})$ for ${\boldsymbol{z}}$ in the window of
${\boldsymbol{z}}_{i}$:

$P_{i}({\boldsymbol{z}})=P_{b,i}({\boldsymbol{z}})e^{\frac{1}{2}\beta\kappa\left|{\boldsymbol{z}}-{\boldsymbol{z}}_{i}\right|^{2}}\left<e^{-\beta\frac{1}{2}\kappa\left|{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}_{i}\right|^{2}}\right>$
(24)

We could use Eq. 24 directly assuming the biased MD trajectory is ergodic, but
we know that regions far from the reference point will be explored very rarely
and thus their free energy would be estimated with large uncertainty. This
means that, although we can use sampling to compute $P_{b,i}$ knowing it
effectively vanishes outside the neighborhood of ${\boldsymbol{z}}_{i}$, we
cannot use sampling to compute
$\left<e^{-\beta\frac{1}{2}\kappa\left|{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}_{i}\right|^{2}}\right>$.

WHAM solves this problem by renormalizing the probabilities in each window
into a single composite probability. Where there is overlap among windows,
WHAM renormalizes such that the statistical variance of the probability is
minimal. That is, it treats the factor
$\left<e^{-\beta\frac{1}{2}\kappa\left|{\boldsymbol{\theta}}({\boldsymbol{x}})-{\boldsymbol{z}}_{i}\right|^{2}}\right>$
as an undetermined constant $C_{i}$ for each window, and solves for specific
values such that the composite unbiased probability $P({\boldsymbol{z}})$ is
continuous across all overlap regions with minimal statistical error. An
alternative to WHAM, termed “umbrella integration”, solves the problem of
renormalization across windows by constructing the composite mean force
Kaestner and Thiel (2005); Kaestner (2011).

The literature on umbrella sampling is vast (by simulation standards), so we
present here a very condensed listing of some of its more recent application
areas with representative citations:

  1. 1.

Small molecule conformational sampling Schaefer et al. (1998); Banavali and
MacKerell (2002); Cruz et al. (2011); Islam et al. (2011);

  2. 2.

Protein-folding Young and Brooks (1996); Sheinerman and Brooks (1998);
Bursulaya and Brooks (1999) and large-scale protein conformational sampling
Rick et al. (1998); Allen et al. (2003); Shams et al. (2012); Yildirim et al.
(2013);

  3. 3.

Protein-protein/peptide-peptide interactions Masunov and Lazaridis (2003);
Tarus et al. (2005); Makowski et al. (2011); Casalini et al. (2011);
Wanasundara et al. (2011); Zhang et al. (2011); Periole et al. (2012);
Vijayaraj et al. (2012); Mahdavi and Kuyucak (2013);

  4. 4.

DNA conformational changes Banavali and Roux (2005) and DNA-DNA interactions
Giudice et al. (2003); Matek et al. (2012); Bagai et al. (2013);

  5. 5.

Binding and association free-energies Czaplewski et al. (2000); Lee and Olson
(2006); Peri et al. (2011); St-Pierre et al. (2011); Rashid and Kuyucak
(2012); Chen and Chung (2012); Wilhelm et al. (2012); Louet et al. (2012);
Zhang et al. (2012); Kessler et al. (2012); Mascarenhas and Kaestner (2013);

  6. 6.

Adsorption on and permeation through lipid bilayers MacCallum and Tieleman
(2006); Tieleman and Marrink (2006); Kyrychenko et al. (2011); Lemkul and
Bevan (2011); Paloncyova et al. (2012); Samanta et al. (2012); Grafmueller et
al. (2013); Cerezo et al. (2013); Tian et al. (2013); Karlsson et al. (2013);

  7. 7.

Adsorption onto inorganic surfaces/interfaces Euston et al. (2011); Doudou et
al. (2012);

  8. 8.

Water ionization Pomes and Roux (1998); Jagoda-Cwiklik et al. (2011);

  9. 9.

Phase transitions Calvo and Mottet (2011); Sharma and Debenedetti (2012);

  10. 10.

Enzymatic mechanisms Ridder et al. (2002); Kaestner et al. (2006); Wang et al.
(2007); Ke et al. (2011); Yan et al. (2011); Mujika et al. (2012); Lonsdale et
al. (2012); Rooklin et al. (2012); Lior-Hoffmann et al. (2012);

  11. 11.

Molecule/ion transport through protein complexes Crouzy et al. (2001); Allen
et al. (2003); Hub and De Groot (2008); Xin et al. (2011); Furini and Domene
(2011); Kim and Allen (2011); Domene and Furini (2012); Zhu and Hummer (2012)
and other macromolecules Zhongjin and Jian (2011); Nalaparaju and Jiang
(2012).

####  2.3.3 Metadynamics

As already mentioned, one of the difficulties of the umbrella sampling method
is the choice and construction of the bias potential. As we already saw with
the relationship among TI, ABF, and TAMD, an adaptive method for building a
bias potential in a running MD simulation may be advantageous. Metadynamics
Laio and Parrinello (2002); Barducci et al. (2011) represents just such a
method.

Metadynamics is rooted in the original idea of “local elevation” Huber et al.
(1994), in which a supplemental bias potential is progressively grown in the
dihedral space of a molecule to prevent it from remaining in one region of
configuration space. However, at variance with metadynamics, local elevation
does not provide any means to reconstruct the unbiased free-energy landscape
and as such it is mostly aimed at fast generation of plausible conformers.

In metadynamics, configurational variables ${\boldsymbol{x}}$ evolve in
response to a biased total potential:

$V({\boldsymbol{x}})=V_{0}({\boldsymbol{x}})+\Delta V({\boldsymbol{x}},t)$
(25)

where $V_{0}$ is the bare interatomic potential and $\Delta
V({\boldsymbol{x}},t)$ is a time-dependent bias potential. The key element of
metadynamics is that the bias is built as a sum of Gaussian functions centered
on the points in CV space already visited:

$\Delta
V\left[{\boldsymbol{\theta}}({\boldsymbol{x}}),t\right]=w\sum_{\begin{array}[]{c}t^{\prime}=\tau_{G},2\tau_{G},\dots\\\
t^{\prime}<t\end{array}}\exp\left(-\frac{\left|{\boldsymbol{\theta}}\left[{\boldsymbol{x}}(t)\right]-{\boldsymbol{\theta}}\left[{\boldsymbol{x}}(t^{\prime})\right]\right|^{2}}{2\delta{\boldsymbol{\theta}}^{2}}\right)$
(26)

Here, $w$ is the height of each Gaussian, $\tau_{G}$ is the size of the time
interval between successive Gaussian depositions, and
$\delta{\boldsymbol{\theta}}$ is the Gaussian width. It has been first
empirically Laio et al. (2005) then analytically Bussi et al. (2006)
demonstrated that in the limit in which the CV evolve according to a Langevin
dynamics, the bias indeed converges to the negative of the free energy, thus
providing an optimal bias to enhance transition events. Multiple simulations
can also be used to allow for a quicker filling of the free-energy landscape
Raiteri et al. (2006).

The difference between the metadynamics estimate of the free energy and the
true free energy can be shown to be related to the diffusion coefficient of
the collective variables and to the rate at which the bias is grown. A
possible way to decrease this error as a simulation progresses is to decrease
the growth rate of the bias. Well-tempered metadynamics Barducci et al. (2008)
used an optimized schedule to decrease the deposition rate of bias by
modulating the Gaussian height:

$w=\omega_{0}\tau_{G}e^{-\frac{\Delta V({\boldsymbol{\theta}},t)}{k_{B}\Delta
T}}$ (27)

Here, $\omega_{0}$ is the initial “deposition rate”, measured Gaussian height
per unit time, and $\Delta T$ is a parameter that controls the  degree to
which the biased trajectory makes excursions away from free-energy minima. It
is possible to show that using well-tempered metadynamics the bias does not
converge to the negative of the free-energy but to a fraction of it, thus
resulting in sampling the CVs at an effectively higher temperature $T+\Delta
T$, where normal metadynamics is recovered for $\Delta T\rightarrow\infty$. We
notice that other deposition schedules can be used aimed, e.g., at maximizing
the number of round-trips in the CV space Singh et al. (2011). Importantly, it
is possible to recover equilibrium Boltzmann statistics of unbiased collective
variables from samples drawn throughout a well-tempered metadynamics
trajectory Bonomi et al. (2009); it does not seem clear that one can do this
from an ABF trajectory. Finally, it is possible to tune the shape of the
Gaussians on the fly using schemes based on the geometric compression of the
phase space or on the variance of the CVs Branduardi et al. (2012).

In the well-tempered ensemble, the parameter $\Delta T$ can be used to tune
the size of the explored region, in a fashion similar to the fictitious
temperature in TAMD. So both TAMD and well-tempered metadynamics can be used
to explore relevant regions of CV space while surmounting relevant free energy
barriers. However, there are important distictions between the two methods.
First, the main source of error in TAMD rests with how well mean-forces are
approximated, and adiabatic separation, realizable only when the auxiliary
variables ${\boldsymbol{z}}$ never move, is the only way to guarantee they are
perfectly accurate. In practical application, TAMD never achieves perfect
adiabatic separation. In contrast, because the deposition rate of decreases as
a well-tempered trajectory progresses, errors related to poor adiabatic
separation are progressively damped. Second, as already mentioned, TAMD alone
cannot report the free energy, but it also is therefore not practically
limited by the dimensionality of CV space; multicomponent gradients are just
as accurately calculated in TAMD as are single-component gradients.
Metadynamics, as a histogram-filling method, must exhaustively sample a finite
region around any point to know the free energy and its gradients are correct,
which can sometimes limit its utility.

Metadynamics is a powerful method whose popularity continues to grow. In
either its original formulation or in more recent variants, metadynamics has
been employed successfully in several fields, some of which we point out below
with some representative examples:

  1. 1.

Chemical reactions Iannuzzi et al. (2003); McGrath et al. (2013);

  2. 2.

Peptide backbone angle sampling Mantz et al. (2009); Leone et al. (2009);
Melis et al. (2009);

  3. 3.

Protein folding Bussi et al. (2006); Gangupomu and Abrams (2010); Berteotti et
al. (2011); Granata et al. (2013);

  4. 4.

Protein aggregation Baftizadeh et al. (2012);

  5. 5.

Molecular docking Gervasio et al. (2005); Soederhjelm et al. (2012);
Limongelli et al. (2013);

  6. 6.

Conformational rearrangement of proteins Sutto and Gervasio (2013 (published
ahead of print);

  7. 7.

Crystal structure prediction Martonak et al. (2006);

  8. 8.

Nucleation and crystal growth Trudu et al. (2006); Stack et al. (2011);

  9. 9.

and proton diffusion Zhang et al. (2012).

###  2.4 Some Comments on Collective Variables

####  2.4.1 The Physical Fidelity of CV-Spaces

Given a potential $V({\boldsymbol{x}})$, any multidimensional CV
${\boldsymbol{\theta}}({\boldsymbol{x}})$ has a mathematically determined free
energy $F({\boldsymbol{z}})$, and in principle the free-energy methods we
describe here (and others) can use and/or compute it. However, this does not
guarantee that $F$ is meaningful, and a poor choice for
${\boldsymbol{\theta}}({\boldsymbol{x}})$ can render the results of even the
most sophisticated free-energy methods useless for understanding the nature of
actual metastable states and the transitions among them. This puts two major
requirements on any CV space:

  1. 1.

Metastable states and transition states must be unambiguously identified as
energetically separate regions in CV space.

  2. 2.

The CV space must not contain hidden barriers.

The first of these may seem obvious: CV’s are chosen to provide a low-
dimensional description of some important process, say a conformational change
or a chemical reaction or a binding event, and one can’t describe a process
without being able to discriminate states. However, it is not always easy to
find CV’s that do this. Even given representative configurations of two
distinct metastable states, standard MD from these two different initial
configurations may sample partially overlapping regions of CV space, making
ambiguous the assignation of an arbitrary configuration to a state. It may be
in this case that the two representative configurations actually belong to the
same state, or that if there are two states, that no matter what CV space is
overlaid, the barrier separating them is so small that, on MD timescales, they
can be considered rapidly exchanging substates of some larger state.

But a third possibility exists: the two MD simulations mentioned above may in
fact represent very different states. The overlap might just be an artifact of
neglecting to include one or more CV’s that are truly necessary to distinguish
those states. If there is a significant free energy barrier along this
neglected variable, an MD simulation will not cross it, yet may still sample
regions in CV space also sampled by an MD simulation launched from the other
side of this hidden barrier. And it is even worse: if TI or umbrella sampling
is used along a pathway in CV space that neglects an important variable, the
free-energy barriers along that pathway might be totally meaningless.

Hidden barriers can be a significant problem in CV-based free-energy
calculations. Generally speaking, one only learns of a hidden barrier after
postulating its existence and testing it with a new calculation. Detecting
them is not straightforward and often involves a good deal of CV space
exploration. Methods such as TAMD and well-tempered metadynamics offer this
capability, but much more work could be done in the automated detection of
hidden barriers and the “right” CV’s (e.g., Das et al. (2006); Perilla and
Woolf (2012); Ceriotti et al. (2011)).

An obvious way of reducing the likelihood of hidden barriers is to use
increase the dimensionality of CV space. TAMD is well-suited to this because
it is a gradient method, but standard metadynamics, because it is a histogram-
filling method, is not. A recent variant of metadynamics termed
“reconnaissance metadynamics” Tribello et al. (2010) does have the capability
of handling high-dimensional CV spaces. In reconnaissance metadynamics, bias
potential kernels are deposited at the CV space points identified as centers
of clusters detected and measured by an on-the-fly clusterization scheme.
These kernels are hyperspherically symmetric but grow as cluster sizes grow
and are able to push a system out of a CV space basin to discover other
basins. As such, reconnaissance metadynamics is an automated way of
identifying free-energy minima in high-dimensional CV spaces. It has been
applied the identification of configurations of small clusters of molecules
Tribello et al. (2011) and identification of protein-ligand binding poses
Soederhjelm et al. (2012).

####  2.4.2 Some Common and Emerging Types of CV’s

There are very few “best practices” codified for choosing CV’s for any given
system. Most CV’s are developed ad hoc based on the processes that
investigators would like to study, for instance, center-of-mass distance
between two molecules for studying binding/unbinding, or torsion angles for
studying conformational changes, or number of contacts for studying order-
disorder transitions. Cartesian coordinates of centers of mass of groups of
atoms are also often used as CV’s, as are functions of these coordinates.

The potential energy $V({\boldsymbol{x}})$ is also an example of a 1-D CV, and
there have been several examples of using it in CV-based enhanced sampling
methods, such as umbrella sampling Bartels and Karplus (1998), metadynamics
Micheletti et al. (2004) well-tempered metadynamics Bonomi and Parrinello
(2010). In a recent work based on steered MD, it has been shown that also
relevant reductions of the potential energy (e.g. the electrostatic
interaction free-energy) can be used as effective CV’s Do et al. (2013). The
basic rationale for enhanced sampling of $V$ is that states with higher
potential energy often correspond to transition states, and one need make no
assumptions about precise physical mechanisms. Key to its successful use as a
CV, as it is for any CV, is a proper accounting for its entropy; i.e., the
classical density-of-states.

Coarse-graining of particle positions onto Eulerian fields was used early on
in enhanced sampling Roitberg and Elber (1991); here, the value of the field
at any Cartesian point is a CV, and the entire field represents a very high-
dimensional CV. This idea has been put to use recently in the “indirect
umbrella sampling” method of Patel et al. Patel et al. (2011) for computing
free energies of solvation, and string method (Sec. 4.2.1) calculations of
lipid bilayer fusion Mueller et al. (2012). In a similar vein, there have been
recent attempts at variables designed to count the recurrency of groups of
atoms positioned according to given templates, such as $\alpha$-helices paired
$\beta$-strands in proteins Pietrucci and Laio (2009).

We finally mention the possibility of building collective variables based on
set of frames which might be available from experimental data or generated by
means of previous MD simulations. Some of these variables are based on the
idea of computing the distances between the present configuration and a set of
precomputed snapshots. These distances, here indicated with $d_{i}$, where $i$
is the index of the snapshot, are then combined to obtain a coarse
representation of the present configuration, which is then used as a CV. As an
example, one might combine the distances as

$s=\frac{\sum_{i}e^{-\lambda d_{i}}i}{\sum_{i}e^{-\lambda d_{i}}}$ (28)

If the parameter $\lambda$ is properly chosen, this function returns a
continuous interpolation between the indexes of the snapshots which are closer
to the present conformation. If the snapshots are disposed along a putative
path connecting two experimental structures, this CV can be used as a path CV
to monitor and bias the progression along the path Branduardi et al. (2007). A
nice feature of path CVs is that it is straighforward to also monitor the
distance from the putative path. The standard way to do it is by looking at
the distance from the closest reference snapshot, which can be approximately
computed with the following continuous function:

$z=-\lambda^{-1}\log\sum_{i}e^{-\lambda d_{i}}$ (29)

This approach, modified to use internal coordinates, was used recently by
Zinovjev et al. to study the aqueous phase reaction of pyruvate to salycilate,
and in the CO bond-breaking/proton transfer in PchB Zinovjev et al. (2012).

A generalization to multidimensional paths (i.e. sheets) can be obtained by
assigning a generic vector $v_{i}$ to each of the precomputed snapshot and
computing its average Spiwok and Králová (2011):

$s=\frac{\sum_{i}e^{-\lambda d_{i}}v_{i}}{\sum_{i}e^{-\lambda d_{i}}}$ (30)

##  3 Tempering Approaches

“Tempering” refers to a class of methods based on increasing the temperature
of an MD system to overcome barriers. Tempering relies on the fact that
according to the Arrhenius law the rate at which activated (barrier-crossing)
events happen is strongly dependent on the temperature. Thus, an annealing
procedure where the system is first heated and then cooled allows one to
produce quickly samples which are largely uncorrelated. The root of all these
ideas indeed lies in the simulated annealing procedure Kirkpatrick et al.
(1983), a well-known method successfully used in many optimization problems.

###  3.1 Simulated tempering

Simulated annealing is a form of Markov-chain Monte Carlo sampling where the
temperature is artificially modified during the simulation. In particular,
sampling is initially done at a temperature high enough that the simulation
can easily overcome high free-energy barriers. Then, the temperature is
decreased as the simulation proceeds, thus smoothly bringing the simulation to
a local energy minimum. In simulated annealing, a critical parameter is the
cooling speed. Indeed, the probability to reach the global minimum grows as
this speed is decreased.

The search for the global minimum can be interpreted in the same way as
sampling an energy landscape at zero temperature. One could thus imagine to
use simulated annealing to generate conformations at, e.g., room temperature
by slowly cooling conformations starting at high temperature. However, the
resulting ensemble will strongly depend on the cooling speed, thus possibly
providing a  biased result. A better approach consists of the  the so-called
simulated tempering methods Marinari and Parisi (1992). Here, a discrete list
of temperatures $T_{i}$, with $i\in 1\dots N$ are chosen _a priori_ ,
typically spanning a range going from the physical temperature of interest to
a temperature which is high enough to overcome all relevant free energy
barriers. (Note that we do not have to stipulate a CV-space in which those
barriers live.) Then, the index $i$, which indicates at which temperature the
system should be simulated, is evolved with time. Two kind of moves are
possible: (a) normal evolution of the system at fixed temperature, which can
be done with a usual Markov Chain Monte Carlo or molecular dynamics and (b)
change of the index $i$ at fixed atomic coordinates. It is easy to show that
the latter can be performed as a Monte Carlo step with acceptance equal to

$\alpha=\min\left(1,\frac{Z_{j}}{Z_{i}}e^{-\frac{U(x)}{k_{B}T_{j}}+\frac{U(x)}{k_{B}T_{i}}}\right)$
(31)

where $i$ and $j$ are the indexes corresponding to the present temperature and
the new one. The weights $Z_{i}$ should be choosen so as to sample
equivalently all the value of $i$. It must be noticed that also within
molecular dynamics simulations only the potential energy usually appears in
the acceptance. This is due to the fact that the velocities are typically
scaled by a factor $\sqrt{\frac{T_{j}}{T_{i}}}$ upon acceptance. This scaling
leads to a cancellation of the contribution to the acceptance coming from the
kinetic energy. Ultimately, this is related to the fact that the ensemble of
velocities is analytically known _a priori_ , such that it is possible to
adapt the velocities to the new temperature instantaneously.

Estimating these  weights $Z_{i}$ is nontrivial and typically requires a
preliminary step. Moreover, if this estimate is poor the system could spend no
time at the physical temperature, thus spoiling the result. Iterative
algorithms for adjusting these weights have been proposed (see e.g. Park and
Pande (2007)). We also observe that since the temperature sets the typical
value of the potential energy, an effect much similar to that of simulated
tempering with adaptive weights can be obtained by performing a metadynamics
simulation using the potential energy as a CV (Sec. 2.4.2).

###  3.2 Parallel tempering

A smart way to alleviate the issue of finding the correct weights is that of
simulating several replicas at the same time Hansmann (1997); Sugita and
Okamoto (1999). Rather that changing the temperature of a single system, the
defining move proposal in parallel tempering consists of a coordinate swap
between two $T$-replicas with acceptance probability

$\alpha=\min\left(1,e^{\left(\frac{1}{k_{B}T_{j}}-\frac{1}{k_{B}T_{i}}\right)\left[U({\boldsymbol{x}}_{i})-U({\boldsymbol{x}}_{j})\right]}\right)$
(32)

This method is the root of a class of techniques collectively known as
“replica exchange” methods, and the latter name is often used as a synonimous
of parallel tempering. Notably, within this framework it is not necessary to
precompute a set of weights. Indeed, the equal time spent by each replica at
each temperature is enforced by the constraint that only pairwise swaps are
allowed. Moreover, parallel tempering has an additional advantage: since the
replicas are weakly coupled and only interact when exchanges are attempted,
they can be simulated on different computers without the need of a very fast
interconnection (provided, of course, that a single replica is small enough to
run on a single node).

The calculation of the acceptance is very cheap as it is based on the
potential energy which is often computed alongside force evaluation. Thus, one
could in theory exploit also a large number of virtual, rejected exchanges so
as to enhance statistical sampling Frenkel (2004); Coluzza and Frenkel (2005).
Since efficiency of parallel tempering simulation can deteriorate if the
stride between subsequent exchanges is too large Sindhikara et al. (2008);
Bussi (2009), a typical recipe is to choose this stride as small as possible,
with the only limitation of avoiding extra costs due to replica
synchronization. One can push this idea further and implement asynchronous
versions of parallel tempering, where overhead related to exchanges is
minimized Gallicchio et al. (2008); Bussi (2009). One should be however aware
that, especially at high exchange rate, artifacts coming from e.g. the use of
wrong thermostating schemes could spoil the results Rosta et al. (2009);
Sindhikara et al. (2010).

Parallel tempering is popular in simulations of protein conformational
sampling Vreede et al. (2005); Zhang and Mu (2012), protein folding Sugita and
Okamoto (1999); Zhou (2003); Garcia and Onuchic (2003); Im et al. (2003); Mei
et al. (2012); Berhanu et al. (2013) and aggregation Kokubo and Okamoto
(2004); Oshaben et al. (2012), due at least in part to the fact that one need
not choose CV’s to use it, and CV’s for describing these processes are not
always straightforward to determine.

###  3.3 Generalized replica exchange

The difference between the replicas is not restricted to be a change in
temperature. Any control parameter can be changed, and even the expression of
the Hamiltonian can be modified Sugita and Okamoto (2000). In the most general
case every replica is simulated at a different temperature (and or pressure)
and a different Hamiltonian, and the acceptance reads

$\alpha=\min\left(1,\frac{e^{-\left(\frac{U_{i}(x_{j})}{k_{B}T_{i}}+\frac{U_{j}(x_{i})}{k_{B}T_{j}}\right)}}{e^{-\left(\frac{U_{i}(x_{i})}{k_{B}T_{i}}+\frac{U_{j}(x_{j})}{k_{B}T_{j}}\right)}}\right)$
(33)

Several recipes for choosing the modified Hamiltonian have been proposed in
the literature Fukunishi et al. (2002); Liu et al. (2005); Affentranger et al.
(2006); Fajer et al. (2008); Xu et al. (2008); Zacharias (2008); Vreede et al.
(2009); Itoh et al. (2010); Meng and Roitberg (2010); Terakawa et al. (2011);
Wang et al. (2011); Zhang and Ma (2012); Bussi (2013). Among these, a notable
idea is that of solute tempering Liu et al. (2005); Wang et al. (2011) which
is used for the simulation of solvated biomolecules. Here, only the
Hamiltonian of the solute is modified. More precisely, one could notice that a
scaling of the Hamiltonian by a factor $\lambda$ is completely equivalent to a
scaling of the temperature by a factor $\lambda^{-1}$. Hamiltonian scaling
however can take advantage of the fact that the total energy of the system is
an extensive property. Thus, one can limit the scaling to the portion of the
system which is considered to be interesting and which has the relevant
bottlenecks. With solute tempering, the solute energy is scaled whereas the
solvent energy is left unchanged. This is equivalent to keeping the solute at
a high effective temperature and the solvent at the physical temperature.
Since in the simulation of solvated molecules most of the atoms belong to the
solvent, this turns in a much smaller modification to the explored ensemble
when compared with parallel tempering. In spite of this, the effect on the
solute resemble much that of increasing the physical temperature.

A sometimes-overlooked subtlety in solute tempering is the choice for the
treatment of solvent-solute interactions. Indeed, whereas solute-solute
interactions are scaled with a factor $\lambda<1$ and solvent-solvent
interactions are not scaled, any intermediate choice (scaling factor between
$\lambda$ and 1) could intuitively make sense for solvent-solute coupling. In
the original formulation, the authors used a factor $(1+\lambda)/2$ for the
solute-solvent interaction. This choice however was later shown to be
suboptimal Huang et al. (2007); Wang et al. (2011), and refined to be
$\sqrt{\lambda}$. This latter choice appears to be more physically sound,
since it allows one to just simulate the biased replicas with a modified
force-field. Indeed, if one scales the charges of the solute by a factor
$\sqrt{\lambda}$, electrostatic interactions are changed by a factor $\lambda$
for solute-solute coupling and $\sqrt{\lambda}$ for solute-solvent coupling.
The same is true for Lennard-Jones terms, albeit in this case it depends on
the specific combination rules used. Notably, the same rules for scaling were
used in a previous work Affentranger et al. (2006). As a final remark, we
point out that solute tempering can be also used in a serial manner _a là_
simulated  tempering, in a simulated solute tempering scheme Denschlag et al.
(2009).

###  3.4 General comments

In general, the advantage of these tempering methods over straighforward
sampling can be rationalized as follows. A simulation is evolved so as to
sample a modified ensemble by e.g. raising temperature or artificially
modifying the Hamiltonian. The change in the ensemble could be drastic, so
that trying to extract canonical averages by reweighting from such a
simulation would be pointless. For this reason, a ladder of intermediate
ensembles is built, interpolating between the physical one (i.e. room
temperature, physical Hamiltonian) and the modified one. Then, transitions
between consecutive steps in this ladder (or, in parallel schemes, coordinate
swaps) are performed using a Monte Carlo scheme. Assuming that the dynamics of
the most modified ensemble is ergodic, independent samples will be generated
every time a new simulation reaches the highest step of the ladder. Thus,
efficiency of these methods is often based on the evaluation of the round trip
time required for a replica to traverse the entire ladder.

Tempering methods are thus relying on the ergodicity of the most modified
ensemble. This assumption is not always correct. A very simple example is
parallel tempering used to accelerate the sampling over an entropic barrier.
Since the height of an entropic barrier grows with the temperature, in this
conditions the barrier in the most modified ensembles are unaffected Zuckerman
and Lyman (2006). Moreover, since a lot of time is spent in sampling states in
non-physical situations (e.g. high temperature), the overall computational
efficiency could even be lower than that of straightforward sampling. Real
applications are often in an intermediate situation, and usefulness of
parallel tempering should be evaluated case by case.

The number of intermediate steps in the ladder can be shown to grow with the
square root of the specific heat of the system in the case of parallel
tempering simulations. No general relationship can be drawn in the case of
Hamiltonian replica exchange, but one can expect approximately that the number
of replicas should be proportional to the square root of the number of degrees
of freedom affected by the modification of the Hamiltonian. Thus, Hamiltonian
replica exchange methods could be much more effective than simple parallel
tempering as they allow the effort to be focused and the number of replicas to
be minimized.

Parallel tempering has the advantage that all the replicas can be analyzed to
obtain meaningful results, e.g., to predict the melting curve of a molecule.
This procedure should be used with caution, especially with empirically
parametrized potentials, which are often tuned to be realistic only at room
temperature. On the other hand, Hamiltonian replica exchange often relies on
unphysically modified ensembles which have no interest but for the fact that
they increase ergodicity.

As a final note, we observe that data obtained at different temperature (or
with modified Hamiltonians) could be combined to enhance statistics at the
physical temperature Chodera et al. (2007). However, the effectiveness of this
data recycling is limited by the fact that high temperature replicas visit
very rarely low energy conformations, thus decreasing the amount of additional
information that can be extracted.

##  4 Combinations and Advanced Approaches

###  4.1 Combination of tempering methods and biased sampling

The algorithms presented in Section 3 and based on tempering are typically
considered to be simpler to apply when compared with those discussed in
Section 2 and based on biasing the sampling of selected collective variables.
Indeed, by avoiding the problem of choosing collective variables which
properly describe the reaction path, most of the burden of setting up a
simulation is removed. However, this comes at a price: considering the
computational cost, tempering methods are extremely expensive. This cost is
related to the fact that they are able to accelerate all degrees of freedom to
the same extent, without an _a priori_ knowledge of the sampling bottlenecks.
In this sense, Hamiltonian replica exchange methods are in an intermediate
situation, since they are typically less expensive than parallel tempering but
allow to embed part of the knowledge of the system in the simulation set up.

Because of the conceptual difference between tempering methods and CV-based
methods, these approaches can be easily and efficiently combined. As an
example, the combination of metadynamics and parallel tempering can be used to
take advantage of the known bottlenecks with biased collective variables at
the same time accelerating the overall sampling with parallel tempering Bussi
et al. (2006). In that work, the free energy landscape for the folding of a
small hairpin was computed by biasing a small number of selected CVs (gyration
radius and the number of hydrogen bonds). These CVs alone are not enough to
describe folding, as can be easily shown by performing a metadynamics
simulation using these CVs. However, the combination with parallel tempering
allowed acceleration of all the degrees of freedom blindly and reversible
folding of the hairpin. This combined approach also improves the results when
compared with parallel tempering alone, since it accelerates exploration of
phase-space. Moreover, since parallel tempering samples the unbiased canonical
distribution, it is very difficult to use it to compute free-energy
differences which are larger than a few $k_{B}T$. The metadynamics bias can be
used to disfavor, e.g., the folded state so as to better estimate the free-
energy difference between the folded and unfolded states.

It is also possible to combine metadynamics with the solute tempering method
so as to decrease the number of required replicas and the computational cost
Camilloni et al. (2008). As an alternative to solute tempering, metadynamics
in the well-tempered ensemble can be effectively used to enhance the
acceptance in parallel tempering simulations and to decrease the number of
necessary replicas Bonomi and Parrinello (2010). This combination of parallel
tempering with well-tempered ensemble can be pushed further and combined with
metadynamics on a few selected degrees of freedom Deighan et al. (2012). As a
final note, bias exchange molecular dynamics Piana and Laio (2007) combines
metadynamics and replica echange in a completely different spirit: every
replica is run using a different CV, thus allowing many CVs to be tried at the
same time. This technique has been succesfully applied to several problems.
For a recent review, we refer the reader to Ref. Baftizadeh et al. (2012).

###  4.2 Some methods based on TAMD

####  4.2.1 String method in collective variables

The string method is generally an approach to find pathways of minimal energy
connecting two points in phase space E et al. (2002). When working in CV’s,
the string method is used to find minimal free-energy paths (MFEP’s)
Maragliano et al. (2006). String method calculations involve multiple
replicas, each representing a point ${\boldsymbol{z}}_{s}$ in CV space at
position $s$ along a discretized string connecting two points of interest
(reactant and product states, say). The forces on each replica’s
${\boldsymbol{z}}_{s}$ are computed and their ${\boldsymbol{z}}_{s}$’s
updated, as in TAMD, with the addition of forces that act to keep the
${\boldsymbol{z}}$’s equidistant along the string (so-called
reparameterization forces):

$\bar{\gamma}\dot{z}_{j}(s,t)=\displaystyle\sum_{k}\biggl{[}\tilde{M}_{jk}(\mathbf{x}(s,t))\kappa[\theta_{k}(\mathbf{x}(s,t))-z_{k}(s,t)]\biggr{]}+\eta_{z}(t)+\lambda(s,t)\displaystyle\frac{\partial{}z_{j}}{\partial{}s}$
(34)

Here, $\tilde{M}_{jk}$ is the metric tensor mapping distances on the manifold
of atomic coordinates to the manifold of CV space, $\eta$ is thermal noise and
$\lambda(s,t)\displaystyle\frac{\partial{}z_{j}}{\partial{}s}$ represents the
reparameterization force tangent to the string that is sufficient to maintain
equidistant images along the string. String method has been used to study
activation of the insulin-receptor kinase Vashisth and Abrams (2012), docking
of insulin to its receptor Vashisth and Abrams (2013), myosin Ovchinnikov et
al. (2011), In these examples, the update of the string coordinates is done at
a lower frequency than the atomic variables in each image.

In contrast, in the on-the-fly variant of string method in CV’s, the friction
on the ${\boldsymbol{z}}_{s}$’s is set high enough to make the effective
averaging of the forces approach the true mean forces, and the
${\boldsymbol{z}}$ updates occur in lockstep with the ${\boldsymbol{x}}$
updates of the MD system Maragliano and Vanden-Eijnden (2007). Just as in
TAMD, the atomic variables obey an equation of motion like Eq. 12 tethering
them to the ${\boldsymbol{z}}_{s}$. Stober and Abrams recently demonstrated an
implementation of on-the-fly string method to study the thermodynamics of the
normal-to-amyloidogenic transition of $\beta$2-microglobulin Stober and Abrams
(2012). Unique in this approach was the construction of a single composite MD
system containing 27 individual $\beta$2 molecules restrained to points on 3
$\times$ 3 $\times$ 3 grid inside a single large solvent box. Zinovjev et al.
used a combination of the on-the-fly string method and of path-collective
variables (see Equations 28 and 29) in a quantum-mechanics/molecular-mechanics
approach to study a methyltransferase reaction Zinovjev et al. (2013).

####  4.2.2 On-the-fly free energy parameterization

Because TAMD provides mean-force estimates as it is exploring CV space, it
stands to reason that those mean forces could be used to compute a free
energy. In contrast, in the single-sweep method Maragliano and Vanden-Eijnden
(2008), the TAMD forces are only used in the CV space exploration phase, not
the free-energy calculation itself. Recently, Abrams and Vanden-Eijnden
proposed a method for using TAMD directly to parameterize a free energy; that
is, to determine the best set of some parameters ${\boldsymbol{\lambda}}$ on
which a free energy of known functional form depends Abrams and Vanden-Eijnden
(2012):

$F({\boldsymbol{z}})=F({\boldsymbol{z}};{\boldsymbol{\lambda}}^{*})$ (35)

The approach, termed “on-the-fly free energy parameterization”, uses forces
from a running TAMD simulation to progressively optimize
${\boldsymbol{\lambda}}$ using a time-averaged gradient error:

$E({\boldsymbol{\lambda}})=\frac{1}{2t}\int_{0}^{t}\left|\nabla_{z}F\left[{\boldsymbol{z}}(s),{\boldsymbol{\lambda}}(t)\right]+\kappa\left[\theta({\boldsymbol{x}}(s))-{\boldsymbol{z}}(s)\right]\right|^{2}ds,$
(36)

If constructed so that  $F$ is linear in
${\boldsymbol{\lambda}}=(\lambda_{1},\lambda_{2},\dots,\lambda_{M})$,
minimization of $E$ can be expressed as a simple linear algebra problem

$\sum_{j}A_{ij}\lambda_{j}=b_{i},\ \ \ i=1,\dots,M$ (37)

and the running TAMD simulation provides progressively better estimates of $A$
and $b$ until the ${\boldsymbol{\lambda}}$ converge. In the cited work, it was
shown that this method is an efficient way to derive potentials of mean force
between particles in coarse-grained molecular simulations as basis-function
expansions. It is currently being investigated as a means to parameterize free
energies associated with conformational changes of proteins.

Chen, Cuendet, and Tuckermann developed a very similar approach that in
addition to parameterizing a free energy using d-AFED-computed gradients uses
a metadynamics-like bias on the potential Chen et al. (2012). These authors
demonstrated efficient reconstruction of the four-dimensional free-energy of
vacuum alanine dipeptide with this approach.

##  5 Concluding Remarks

In this review, we have summarized some of the current and emerging enhanced
sampling methods that sit atop MD simulation. These have been broadly
classified as methods that use collective variable biasing and methods that
use tempering. CV biasing is a much more prevalent approach than tempering,
due partially to the fact that it is perceived to be cheaper, since tempering
simulations are really only useful for enhanced sampling of configuration
space when run in parallel. CV-biasing also reflects the desire to rein in the
complexity of all-atom simulations by projecting configurations into a much
lower dimensional space. (Parallel tempering can be thought of as increasing
the dimensionality of the system by a factor equal to the number of simulated
replicas.) But the drawback of all CV-biasing approaches is the risk that the
chosen CV space does not provide the most faithful representation of the true
spectrum of metastable subensembles and the barriers that separate them.
Guaranteeing that sampling of CV space is not stymied by hidden barriers must
be of paramount concern in the continued evolution of such methods. For this
reason, methods that specifically allow broad exploration of CV space, like
TAMD (which can handle large numbers of CV’s) and well-tempered metadynamics
will continue to be valuable. So too will parallel tempering because its broad
sampling of configuration space can be used to inform the choice of better
CV’s. Accelerating development of combined CV-tempering methods bodes well for
enhanced sampling generally.

Although some of these methods involve time-varying forces (ABF, TAMD, and
metadynamics), all methods we’ve discussed have the underlying rationale of
the equilibrium ensemble. TI uses the constrained ensemble, ABF and
metadynamics ideally converge to an ensemble in which a bias erases free-
energy variations, and TAMD samples an attenuated/mollified equilibrium
ensemble. There is an entirely separate class of methods that inherently rely
on non-equilibrium thermodynamics. We have not discussed at all the several
free-energy methods based on non-equilibrium MD simulations; we refer
interested readers to the article by Christoph Dellago in this issue.

Finally, we have also not really touched on any of the practical issues of
implementing and using these methods in conjunction with modern MD packages
(e.g., NAMD Phillips et al. (2005), LAMMPS Plimpton (1995), Gromacs Hess et
al. (2008), Amber Case et al. (2005), and CHARMM Brooks et al. (1983), to name
a few). At least two packages (NAMD and CHARMM) have native support for
collective variable biasing, and NAMD in particular offers both native ABF and
a TcL-based interface which has been used to implement TAMD Abrams and Vanden-
Eijnden (2010). The native collective variable module for NAMD has been
recently ported to LAMMPS Fiorin et al. (2013). Gromacs offers native support
for parallel tempering. Generally speaking, however, modifying MD codes to
handle CV-biasing and multiple replicas is not straightforward, since one
would like access to the data structures that store coordinates and forces. A
major help in this regard is the PLUMED package Bonomi et al. (2009); Tribello
et al. (2013), which patches a variety of MD codes to enable users to use many
of the techniques discussed here.

##  6 Abbreviations

  * •

ABF: adaptive-biasing force

  * •

AFED: adiabatic free-energy dynamics

  * •

CV: collective-variable

  * •

MD: molecular dynamics

  * •

MFEP: minimum free-energy path

  * •

TAMD: temperature-accelerated molecular dynamics

  * •

TI: thermodynamic integration

  * •

WHAM: weighted-histogram analysis method

## Acknowledgments

CFA would like to acknowledge support of NSF (DMR-1207389) and NIH
(1R01GM100472). GB would like to acknowledge the European Research Council
(Starting Grant S-RNA-S, no. 306662) for financial support. Both authors would
like to acknowledge NSF support of a recent Pan-American Advanced Studies
Institute Workshop “Molecular-based Multiscale Modeling and Simulation”
(OISE-1124480; PI: W. J. Pfaednter, U. Washington) held in Montevideo,
Uruguay, Sept. 12-15, 2012, where the authors met and began discussions that
influenced the content of this review.

## References

  * Kollman (1993) Kollman, P.  Free-energy calculations - applications to chemical and biochemical phenomena.  Chem. Rev. 1993, 93, 2395–2417.
  * Trzesniak et al. (2007) Trzesniak, D.; Kunz, A.P.E.; van Gunsteren, W.F.  A comparison of methods to compute the potential of mean force.  ChemPhysChem 2007, 8, 162–169.
  * Vanden-Eijnden (2009) Vanden-Eijnden, E.  Some Recent Techniques for Free Energy Calculations.  J. Comput. Chem. 2009, 30, 1737–1747.
  * Dellago and Bolhuis (2009) Dellago, C.; Bolhuis, P.G.  Transition path sampling and other advanced simulation techniques for rare events. In Advanced computer simulation approaches for soft matter sciences III; Springer, 2009; pp. 167–233.
  * Christ et al. (2010) Christ, C.D.; Mark, A.E.; van Gunsteren, W.F.  Basic Ingredients of Free Energy Calculations: A Review.  J. Comput. Chem. 2010, 31, 1569–1582.
  * Darve et al. (2008) Darve, E.; Rodriguez-Gomez, D.; Pohorille, A.  Adaptive biasing force method for scalar and vector free energy calculations.  J. Chem. Phys. 2008, 128.
  * Maragliano and Vanden-Eijnden (2006) 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, 168–175.
  * Laio and Parrinello (2002) Laio, A.; Parrinello, M.  Escaping free-energy minima.  Proc. Natl. Acad. Sci. USA 2002, 99, 12562–12566.
  * Schlitter et al. (1993) Schlitter, J.; Engels, M.; Kruger, P.; Jacoby, E.; Wollmer, A.  Targeted molecular-dynamics simulation of conformational change - application to the T[–]T transition in insulin.  Mol. Sim. 1993, 10, 291–&.
  * Grubmüller et al. (1996) Grubmüller, H.; Heymann, B.; Tavan, P.  Ligand binding: Molecular mechanics calculation of the streptavidin biotin rupture force.  Science 1996, 271, 997–999.
  * Kirkwood (1935) Kirkwood, J.G.  Statistical Mechanics of Fluid Mixtures.  J. Chem. Phys. 1935, 3, 300–313.
  * Carter et al. (1989) Carter, E.; Ciccotti, G.; Hynes, J.T.; Kapral, R.  Constrained reaction coordinate dynamics for the simulation of rare events.  Chem. Phys. Lett. 1989, 156, 472 – 477.
  * Sprik and Ciccotti (1998) Sprik, M.; Ciccotti, G.  Free energy from constrained molecular dynamics.  J. Chem. Phys. 1998, 109, 7737–7744.
  * Ciccotti and Ferrario (2004) Ciccotti, G.; Ferrario, M.  Blue moon approach to rare events.  Mol. Sim. 2004, 30, 787–793.
  * Ciccotti et al. (2005) Ciccotti, G.; Kapral, R.; Vanden-Eijnden, E.  Blue moon sampling, vectorial reaction coordinates, and unbiased constrained dynamics.  ChemPhysChem 2005, 6, 1809–1814.
  * den Otter and Briels (1998) den Otter, W.; Briels, W.  The calculation of free-energy differences by constrained molecular-dynamics simulations.  J. Chem. Phys. 1998, 109, 4139–4146.
  * Schlitter and Klähn (2003) Schlitter, J.; Klähn, M.  A new concise expression for the free energy of a reaction coordinate.  J. Chem. Phys. 2003, 118, 2057–2060.
  * Depaepe et al. (1993) Depaepe, J.; Ryckaert, J.; Paci, E.; Ciccotti, G.  Sampling of molecular-conformations by molecular-dynamics techniques.  Mol. Phys. 1993, 79, 515–522.
  * Zhao et al. (2008) Zhao, X.; Rignall, T.R.; McCabe, C.; Adney, W.S.; Himmel, M.E.  Molecular simulation evidence for processive motion of Trichoderma reesei Cel7A during cellulose depolymerization.  Chem. Phys. Lett. 2008, 460, 284–288.
  * Kim et al. (2009) Kim, H.; Goddard, III, W.A.; Jang, S.S.; Dichtel, W.R.; Heath, J.R.; Stoddart, J.F.  Free Energy Barrier for Molecular Motions in Bistable [2]Rotaxane Molecular Electronic Devices.  J. Phys. Chem. A 2009, 113, 2136–2143.
  * Hytha et al. (2001) Hytha, M.; Stich, I.; Gale, J.; Terakura, K.; Payne, M.  Thermodynamics of catalytic formation of dimethyl ether from methanol in acidic zeolites.  Chem.-A Eur. J. 2001, 7, 2521–2527.
  * Fois et al. (2004) Fois, E.; Gamba, A.; Spano, E.  Competition between water and hydrogen peroxide at Ti center in Titanium zeolites. An ab initio study.  J. Phys. Chem. B 2004, 108, 9557–9560.
  * Ivanov and Klein (2005) Ivanov, I.; Klein, M.  Dynamical flexibility and proton transfer in the arginase active site probed by ab initio molecular dynamics.  J. Am. Chem. Soc. 2005, 127, 4010–4020.
  * Stubbs and Marx (2005) Stubbs, J.; Marx, D.  Aspects of glycosidic bond formation in aqueous solution: Chemical bonding and the role of water.  Chem.-A Eur. J. 2005, 11, 2651–2659.
  * Trinh et al. (2009) Trinh, T.T.; Jansen, A.P.J.; van Santen, R.A.; Meijer, E.J.  The role of water in silicate oligomerization reaction.  Phys. Chem. Chem. Phys. 2009, 11, 5092–5099.
  * Liu et al. (2010) Liu, P.; Cai, W.; Chipot, C.; Shao, X.  Thermodynamic Insights into the Dynamic Switching of a Cyclodextrin in a Bistable Molecular Shuttle.  J. Phys. Chem. Lett. 2010, 1, 1776–1780.
  * Bucko and Hafner (2010) Bucko, T.; Hafner, J.  Entropy effects in hydrocarbon conversion reactions: free-energy integrations and transition-path sampling.  J. Phys.-Cond. Mat. 2010, 22.
  * Paci and Marchi (1994) Paci, E.; Marchi, M.  Membrane crossing by a polar molecule - a molecular-dynamics simulation.  Mol. Sim. 1994, 14, 1–10.
  * Sa et al. (2006) Sa, R.; Zhu, W.; Shen, J.; Gong, Z.; Cheng, J.; Chen, K.; Jiang, H.  How does ammonium dynamically interact with benzene in aqueous media? A first principle study using the Car-Parrinello molecular dynamics method.  J. Phys. Chem. B 2006, 110, 5094–5098.
  * Mugnai et al. (2007) Mugnai, M.; Cardini, G.; Schettino, V.; Nielsen, C.J.  Ab initio molecular dynamics study of aqueous formaldehyde and methanediol.  Mol. Phys. 2007, 105, 2203–2210.
  * Chunsrivirot and Trout (2011) Chunsrivirot, S.; Trout, B.L.  Free energy of binding of a small molecule to an amorphous polymer in a solvent.  Langmuir 2011, 27, 6910–6919.
  * Sato et al. (2012) Sato, M.; Yamataka, H.; Komeiji, Y.; Mochizuki, Y.  FMO-MD Simulations on the Hydration of Formaldehyde in Water Solution with Constraint Dynamics.  Chem.-A Eur. J. 2012, 18, 9714–9721.
  * Sergi et al. (2002) Sergi, A.; Ciccotti, G.; Falconi, M.; Desideri, A.; Ferrario, M.  Effective binding force calculation in a dimeric protein by molecular dynamics simulation.  J. Chem. Phys. 2002, 116, 6329–6338.
  * Maragliano et al. (2004) Maragliano, L.; Ferrario, M.; Ciccotti, G.  Effective binding force calculation in dimeric proteins.  Mol. Sim. 2004, 30, 807–816.
  * Darve and Pohorille (2001) Darve, E.; Pohorille, A.  Calculating free energies using average force.  J. Chem. Phys. 2001, 115, 9169–9183.
  * Fogolari et al. (2011) Fogolari, F.; Corazza, A.; Varini, N.; Rotter, M.; Gumral, D.; Codutti, L.; Rennella, E.; Viglino, P.; Bellotti, V.; Esposito, G.  Molecular dynamics simulation of beta(2)-microglobulin in denaturing and stabilizing conditions.  Proteins 2011, 79, 986–1001.
  * Faller et al. (2013) Faller, C.E.; Reilly, K.A.; Hills, Jr., R.D.; Guvench, O.  Peptide Backbone Sampling Convergence with the Adaptive Biasing Force Algorithm.  J. Phys. Chem. B 2013, 117, 518–526.
  * Wei and Pohorile (2011) Wei, C.; Pohorile, A.  Permeation of Nucleosides through Lipid Bilayers.  J. Phys. Chem. B 2011, 115, 3681–3688.
  * Vivcharuk and Kaznessis (2011) Vivcharuk, V.; Kaznessis, Y.N.  Thermodynamic Analysis of Protegrin-1 Insertion and Permeation through a Lipid Bilayer.  J. Phys. Chem. B 2011, 115, 14704–14712.
  * Kraszewski et al. (2011) Kraszewski, S.; Tarek, M.; Ramseyer, C.  Uptake and Translocation Mechanisms of Cationic Amino Derivatives Functionalized on Pristine C-60 by Lipid Membranes: A Molecular Dynamics Simulation Study.  ACS Nano 2011, 5, 8571–8578.
  * Kraszewski et al. (2012) Kraszewski, S.; Bianco, A.; Tarek, M.; Ramseyer, C.  Insertion of Short Amino-Functionalized Single-Walled Carbon Nanotubes into Phospholipid Bilayer Occurs by Passive Diffusion.  PLOS ONE 2012, 7.
  * Liu et al. (2010) Liu, X.; Lu, X.; Meijer, E.J.; Wang, R.; Zhou, H.  Acid dissociation mechanisms of Si(OH)(4) and Al(H2O)(6)(3+) in aqueous solution.  Geochim. Cosmochim. Acta 2010, 74, 510–516.
  * Caballero et al. (2013) Caballero, J.; Poblete, H.; Navarro, C.; Alzate-Morales, J.H.  Association of nicotinic acid with a poly(amidoamine) dendrimer studied by molecular dynamics simulations.  J. Molec. Graph. Mod. 2013, 39, 71–78.
  * Wilson et al. (2011) Wilson, M.A.; Wei, C.; Bjelkmar, P.; Wallace, B.A.; Pohorille, A.  Molecular Dynamics Simulation of the Antiamoebin Ion Channel: Linking Structure and Conductance.  Biophys. J. 2011, 100, 2394–2402.
  * Cheng and Coalson (2012) Cheng, M.H.; Coalson, R.D.  Molecular Dynamics Investigation of Cl- and Water Transport through a Eukaryotic CLC Transporter.  Biophys. J. 2012, 102, 1363–1371.
  * Wang et al. (2012) Wang, S.; Orabi, E.A.; Baday, S.; Berneche, S.; Lamoureux, G.  Ammonium Transporters Achieve Charge Transfer by Fragmenting Their Substrate.  J. Am. Chem. Soc. 2012, 134, 10419–10427.
  * Tillman et al. (2013) Tillman, T.; Cheng, M.H.; Chen, Q.; Tang, P.; Xu, Y.  Reversal of ion-charge selectivity renders the pentameric ligand-gated ion channel GLIC insensitive to anaesthetics.  Biochem. J. 2013, 449, 61–68.
  * Akhshi et al. (2012) Akhshi, P.; Acton, G.; Wu, G.  Molecular Dynamics Simulations to Provide New Insights into the Asymmetrical Ammonium Ion Movement Inside of the [d(G(3)T(4)G(4))](2) G-Quadruplex DNA Structure.  J. Phys. Chem. B 2012, 116, 9363–9370.
  * Kamath et al. (2012) Kamath, G.; Bhatnagar, N.; Baker, G.A.; Baker, S.N.; Potoff, J.J.  Computational prediction of ionic liquid 1-octanol/water partition coefficients.  Phys. Chem. Chem. Phys. 2012, 14, 4339–4342.
  * Bhatnagar et al. (2012) Bhatnagar, N.; Kamath, G.; Chelst, I.; Potoff, J.J.  Direct calculation of 1-octanol-water partition coefficients from adaptive biasing force molecular dynamics simulations.  J. Chem. Phys. 2012, 137.
  * Wereszczynski and McCammon (2012) Wereszczynski, J.; McCammon, J.A.  Nucleotide-dependent mechanism of Get3 as elucidated from free energy calculations.  Proc. Natl. Acad. Sci. USA 2012, 109, 7759–7764.
  * Jana and Sengupta (2012) Jana, A.K.; Sengupta, N.  Adsorption Mechanism and Collapse Propensities of the Full-Length, Monomeric A beta(1-42) on the Surface of a Single-Walled Carbon Nanotube: A Molecular Dynamics Simulation Study.  Biophys. J. 2012, 102, 1889–1896.
  * Uddin et al. (2010) Uddin, N.M.; Capaldi, F.; Farouk, B.  Molecular Dynamics Simulations of Carbon Nanotube Interactions in Water/Surfactant Systems.  J. Eng. Mater.-T. ASME 2010, 132.
  * Kaestner (2009) Kaestner, J.  Umbrella integration in two or more reaction coordinates.  J. Chem. Phys. 2009, 131, 034109.
  * Rosso and Tuckerman (2002) Rosso, L.; Tuckerman, M.  An adiabatic molecular dynamics method for the calculation of free energy profiles.  Mol. Sim. 2002, 28, 91–112.
  * Abrams and Tuckerman (2008) Abrams, J.B.; Tuckerman, M.E.  Efficient and Direct Generation of Multidimensional Free Energy Surfaces via Adiabatic Dynamics without Coordinate Transformations.  J. Phys. Chem. B 2008, 112, 15742–15757.
  * Iannuzzi et al. (2003) Iannuzzi, M.; Laio, A.; Parrinello, M.  Efficient exploration of reactive potential energy surfaces using Car-Parrinello molecular dynamics.  Phys. Rev. Lett. 2003, 90.
  * Abrams and Vanden-Eijnden (2010) Abrams, C.; Vanden-Eijnden, E.  Large-scale conformational sampling of proteins using temperature-accelerated molecular dynamics.  Proc. Natl. Acad. Sci. USA 2010, 107, 4961–4966.
  * Vashisth et al. (2012) Vashisth, H.; Skiniotis, G.; Brooks, III, C.L.  Using Enhanced Sampling and Structural Restraints to Refine Atomic Structures into Low-Resolution Electron Microscopy Maps.  Structure 2012, 20, 1453–1462.
  * Vashisth et al. (2013) Vashisth, H.; Skiniotis, G.; Brooks, III, C.L.  Enhanced Sampling and Overfitting Analyses in Structural Refinement of Nucleic Acids into Electron Microscopy Maps.  J. Phys. Chem. B 2013, 117, 3738–3746.
  * Vashisth and Brooks (2012) Vashisth, H.; Brooks, III, C.L.  Conformational Sampling of Maltose-Transporter Components in Cartesian Collective Variables Is Governed by the Low-Frequency Normal Modes.  J. Phys. Chem. Lett. 2012, 3, 3379–3384.
  * Hu et al. (2012) Hu, Y.; Hong, W.; Shi, Y.; Liu, H.  Temperature-Accelerated Sampling and Amplified Collective Motion with Adiabatic Reweighting to Obtain Canonical Distributions and Ensemble Averages.  J. Chem. Theor. Comput. 2012, 8, 3777–3792.
  * Vashisth and Abrams (2012) Vashisth, H.; Abrams, C.F.  DFG-flip in the insulin receptor kinease is facilitated by a helical intermediate state of the activation loop.  Biophys. J. 2012, 102, 1979–1987.
  * Maragliano et al. (2010) Maragliano, L.; Cottone, G.; Ciccotti, G.; Vanden-Eijnden, E.  Mapping the Network of Pathways of CO Diffusion in Myoglobin.  J. Am. Chem. Soc. 2010, 132, 1010–1017.
  * Lapelosa and Abrams (2013) Lapelosa, M.; Abrams, C.F.  A Computational Study of Water and CO Migration Sites and Channels Inside Myoglobin.  J. Chem. Theory Comput. 2013, 9, 1265–1271.
  * Geslin et al. (2013) Geslin, P.A.; Ciccotti, G.; Meloni, S.  An observable for vacancy characterization and diffusion in crystals.  J. Chem. Phys. 2013, 138.
  * Lucid et al. (2013) Lucid, J.; Meloni, S.; MacKernan, D.; Spohr, E.; Ciccotti, G.  Probing the Structures of Hydrated Nafion in Different Morphologies Using Temperature-Accelerated Molecular Dynamics Simulations.  J. Phys. Chem. C 2013, 117, 774–782.
  * Maragliano and Vanden-Eijnden (2008) Maragliano, L.; Vanden-Eijnden, E.  Single-sweep methods for free energy calculations.  J. Chem. Phys. 2008, 128.
  * Torrie and Valleau (1977) Torrie, G.M.; Valleau, J.P.  Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling.  J. Comput. Phys. 1977, 23, 187–199.
  * Kumar et al. (1992) Kumar, S.; Rosenberg, J.M.; Bouzida, D.; Swendsen, R.H.; Kollman, P.A.  The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method.  J. Comput. Chem. 1992, 13, 1011–1021.
  * Kaestner and Thiel (2005) Kaestner, J.; Thiel, W.  Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: “umbrella integration”.  J. Chem. Phys. 2005, 123, 144104.
  * Kaestner (2011) Kaestner, J.  Umbrella sampling.  Wires. Comput. Mol. Sci. 2011, 1, 932–942.
  * Schaefer et al. (1998) Schaefer, M.; Bartels, C.; Karplus, M.  Solution conformations and thermodynamics of structured peptides: Molecular dynamics simulation with an implicit solvation model.  J. Mol. Biol. 1998, 284, 835–848.
  * Banavali and MacKerell (2002) Banavali, N.; MacKerell, A.  Free energy and structural pathways of base flipping in a DNA GCGC containing sequence.  J. Mol. Biol. 2002, 319, 141–160.
  * Cruz et al. (2011) Cruz, V.; Ramos, J.; Martinez-Salazar, J.  Water-Mediated Conformations of the Alanine Dipeptide as Revealed by Distributed Umbrella Sampling Simulations, Quantum Mechanics Based Calculations, and Experimental Data.  J. Phys. Chem. B 2011, 115, 4880–4886.
  * Islam et al. (2011) Islam, S.M.; Richards, M.R.; Taha, H.A.; Byrns, S.C.; Lowary, T.L.; Roy, P.N.  Conformational Analysis of Oligoarabinofuranosides: Overcoming Torsional Barriers with Umbrella Sampling.  J. Chem. Theory Comput. 2011, 7, 2989–3000.
  * Young and Brooks (1996) Young, W.; Brooks, C.  A microscopic view of helix propagation: N and C-terminal helix growth in alanine helices.  J. Mol. Biol. 1996, 259, 560–572.
  * Sheinerman and Brooks (1998) Sheinerman, F.; Brooks, C.  Calculations on folding of segment B1 of streptococcal protein G.  J. Mol. Biol. 1998, 278, 439–456.
  * Bursulaya and Brooks (1999) Bursulaya, B.; Brooks, C.  Folding free energy surface of a three-stranded beta-sheet protein.  J. Am. Chem. Soc. 1999, 121, 9947–9951.
  * Rick et al. (1998) Rick, S.; Erickson, J.; Burt, S.  Reaction path and free energy calculations of the transition between alternate conformations of HIV-1 protease.  Proteins 1998, 32, 7–16.
  * Allen et al. (2003) Allen, T.; Andersen, O.; Roux, B.  Structure of gramicidin A in a lipid bilayer environment determined using molecular dynamics simulations and solid-state NMR data.  J. Am. Chem. Soc. 2003, 125, 9868–9877.
  * Shams et al. (2012) Shams, H.; Golji, J.; Mofrad, M.R.K.  A Molecular Trajectory of alpha-Actinin Activation.  Biophys. J. 2012, 103, 2050–2059.
  * Yildirim et al. (2013) Yildirim, I.; Park, H.; Disney, M.D.; Schatz, G.C.  A Dynamic Structural Model of Expanded RNA CAG Repeats: A Refined X-ray Structure and Computational Investigations Using Molecular Dynamics and Umbrella Sampling Simulations.  J. Am. Chem. Soc. 2013, 135, 3528–3538.
  * Masunov and Lazaridis (2003) Masunov, A.; Lazaridis, T.  Potentials of mean force between ionizable amino acid side chains in water.  J. Am. Chem. Soc. 2003, 125, 1722–1730.
  * Tarus et al. (2005) Tarus, B.; Straub, J.; Thirumalai, D.  Probing the initial stage of aggregation of the A beta(10-35)-protein: Assessing the propensity for peptide dimerization.  J. Mol. Biol. 2005, 345, 1141–1156.
  * Makowski et al. (2011) Makowski, M.; Liwo, A.; Sobolewski, E.; Scheraga, H.A.  Simple Physics-Based Analytical Formulas for the Potentials of Mean Force of the Interaction of Amino-Acid Side Chains in Water. V. Like-Charged Side Chains.  J. Phys. Chem. B 2011, 115, 6119–6129.
  * Casalini et al. (2011) Casalini, T.; Salvalaglio, M.; Perale, G.; Masi, M.; Cavallotti, C.  Diffusion and Aggregation of Sodium Fluorescein in Aqueous Solutions.  J. Phys. Chem. B 2011, 115, 12896–12904.
  * Wanasundara et al. (2011) Wanasundara, S.N.; Krishnamurthy, V.; Chung, S.H.  Free Energy Calculations of Gramicidin Dimer Dissociation.  J. Phys. Chem. B 2011, 115, 13765–13770.
  * Zhang et al. (2011) Zhang, B.W.; Brunetti, L.; Brooks, III, C.L.  Probing pH-Dependent Dissociation of HdeA Dimers.  J. Am. Chem. Soc. 2011, 133, 19393–19398.
  * Periole et al. (2012) Periole, X.; Knepp, A.M.; Sakmar, T.P.; Marrink, S.J.; Huber, T.  Structural Determinants of the Supramolecular Organization of G Protein-Coupled Receptors in Bilayers.  J. Am. Chem. Soc. 2012, 134, 10959–10965.
  * Vijayaraj et al. (2012) Vijayaraj, R.; Van Damme, S.; Bultinck, P.; Subramanian, V.  Molecular Dynamics and Umbrella Sampling Study of Stabilizing Factors in Cyclic Peptide-Based Nanotubes.  J. Phys. Chem. B 2012, 116, 9922–9933.
  * Mahdavi and Kuyucak (2013) Mahdavi, S.; Kuyucak, S.  Why the Drosophila Shaker K+ Channel Is Not a Good Model for Ligand Binding to Voltage-Gated Kv1 Channels.  Biochemistry 2013, 52, 1631–1640.
  * Banavali and Roux (2005) Banavali, N.; Roux, B.  Free energy landscape of A-DNA to B-DNA conversion in aqueous solution.  J. Am. Chem. Soc. 2005, 127, 6866–6876.
  * Giudice et al. (2003) Giudice, E.; Varnai, P.; Lavery, R.  Base pair opening within B-DNA: free energy pathways for GC and AT pairs from umbrella sampling simulations.  Nucl. Acids Res. 2003, 31, 1434–1443.
  * Matek et al. (2012) Matek, C.; Ouldridge, T.E.; Levy, A.; Doye, J.P.K.; Louis, A.A.  DNA Cruciform Arms Nucleate through a Correlated but Asynchronous Cooperative Mechanism.  J. Phys. Chem. B 2012, 116, 11616–11625.
  * Bagai et al. (2013) Bagai, S.; Sun, C.; Tang, T.  Potential of Mean Force of Polyethylenimine-Mediated DNA Attraction.  J. Phys. Chem. B 2013, 117, 49–56.
  * Czaplewski et al. (2000) Czaplewski, C.; Rodziewicz-Motowidlo, S.; Liwo, A.; Ripoll, D.; Wawak, R.; Scheraga, H.  Molecular simulation study of cooperativity in hydrophobic association.  Protein Sci. 2000, 9, 1235–1245.
  * Lee and Olson (2006) Lee, M.; Olson, M.  Calculation of absolute protein-ligand binding affinity using path and endpoint approaches.  Biophys. J. 2006, 90, 864–877.
  * Peri et al. (2011) Peri, S.; Karim, M.N.; Khare, R.  Potential of mean force for separation of the repeating units in cellulose and hemicellulose.  Carbohyd. Res. 2011, 346, 867–871.
  * St-Pierre et al. (2011) St-Pierre, J.F.; Karttunen, M.; Mousseau, N.; Rog, T.; Bunker, A.  Use of Umbrella Sampling to Calculate the Entrance/Exit Pathway for Z-Pro-Prolinal Inhibitor in Prolyl Oligopeptidase.  J. Chem. Theory Comput. 2011, 7, 1583–1594.
  * Rashid and Kuyucak (2012) Rashid, M.H.; Kuyucak, S.  Affinity and Selectivity of ShK Toxin for the Kv1 Potassium Channels from Free Energy Simulations.  J. Phys. Chem. B 2012, 116, 4812–4822.
  * Chen and Chung (2012) Chen, R.; Chung, S.H.  Conserved Functional Surface of Antimammalian Scorpion beta-Toxins.  J. Phys. Chem. B 2012, 116, 4796–4800.
  * Wilhelm et al. (2012) Wilhelm, M.; Mukherjee, A.; Bouvier, B.; Zakrzewska, K.; Hynes, J.T.; Lavery, R.  Multistep Drug Intercalation: Molecular Dynamics and Free Energy Studies of the Binding of Daunomycin to DNA.  J. Am. Chem. Soc. 2012, 134, 8588–8596.
  * Louet et al. (2012) Louet, M.; Martinez, J.; Floquet, N.  GDP Release Preferentially Occurs on the Phosphate Side in Heterotrimeric G-proteins.  PLOS Comput. Biol. 2012, 8.
  * Zhang et al. (2012) Zhang, H.; Tan, T.; Feng, W.; van der Spoel, D.  Molecular Recognition in Different Environments: beta-Cyclodextrin Dimer Formation in Organic Solvents.  J. Phys. Chem. B 2012, 116, 12684–12693.
  * Kessler et al. (2012) Kessler, J.; Jakubek, M.; Dolensky, B.; Bour, P.  Binding energies of five molecular pincers calculated by explicit and implicit solvent models.  J. Comput. Chem. 2012, 33, 2310–2317.
  * Mascarenhas and Kaestner (2013) Mascarenhas, N.M.; Kaestner, J.  How maltose influences structural changes to bind to maltose-binding protein: Results from umbrella sampling simulation.  Proteins 2013, 81, 185–198.
  * MacCallum and Tieleman (2006) MacCallum, J.; Tieleman, D.  Computer simulation of the distribution of hexane in a lipid bilayer: Spatially resolved free energy, entropy, and enthalpy profiles.  J. Am. Chem. Soc. 2006, 128, 125–130.
  * Tieleman and Marrink (2006) Tieleman, D.P.; Marrink, S.J.  Lipids out of equilibrium: Energetics of desorption and pore mediated flip-flop.  J. Am. Chem. Soc. 2006, 128, 12462–12467.
  * Kyrychenko et al. (2011) Kyrychenko, A.; Sevriukov, I.Y.; Syzova, Z.A.; Ladokhin, A.S.; Doroshenko, A.O.  Partitioning of 2,6-Bis(1H-Benzimidazol-2-yl)pyridine fluorophore into a phospholipid bilayer: Complementary use of fluorescence quenching studies and molecular dynamics simulations.  Biophys. Chem. 2011, 154, 8–17.
  * Lemkul and Bevan (2011) Lemkul, J.A.; Bevan, D.R.  Characterization of Interactions between PilA from Pseudomonas aeruginosa Strain K and a Model Membrane.  J. Phys. Chem. B 2011, 115, 8004–8008.
  * Paloncyova et al. (2012) Paloncyova, M.; Berka, K.; Otyepka, M.  Convergence of Free Energy Profile of Coumarin in Lipid Bilayer.  J. Chem. Theory Comput. 2012, 8, 1200–1211.
  * Samanta et al. (2012) Samanta, S.; Hezaveh, S.; Milano, G.; Roccatano, D.  Diffusion of 1,2-Dimethoxyethane and 1,2-Dimethoxypropane through Phosphatidycholine Bilayers: A Molecular Dynamics Study.  J. Phys. Chem. B 2012, 116, 5141–5151.
  * Grafmueller et al. (2013) Grafmueller, A.; Lipowsky, R.; Knecht, V.  Effect of tension and curvature on the chemical potential of lipids in lipid aggregates.  Phys. Chem. Chem. Phys. 2013, 15, 876–881.
  * Cerezo et al. (2013) Cerezo, J.; Zuniga, J.; Bastida, A.; Requena, A.; Ceron-Carrasco, J.P.  Conformational changes of beta-carotene and zeaxanthin immersed in a model membrane through atomistic molecular dynamics simulations.  Phys. Chem. Chem. Phys. 2013, 15, 6527–6538.
  * Tian et al. (2013) Tian, J.; Sethi, A.; Swanson, B.I.; Goldstein, B.; Gnanakaran, S.  Taste of Sugar at the Membrane: Thermodynamics and Kinetics of the Interaction of a Disaccharide with Lipid Bilayers.  Biophys. J. 2013, 104, 622–632.
  * Karlsson et al. (2013) Karlsson, B.C.G.; Olsson, G.D.; Friedman, R.; Rosengren, A.M.; Henschel, H.; Nicholls, I.A.  How Warfarin’s Structural Diversity Influences Its Phospholipid Bilayer Membrane Permeation.  J. Phys. Chem. B 2013, 117, 2384–2395.
  * Euston et al. (2011) Euston, S.R.; Bellstedt, U.; Schillbach, K.; Hughes, P.S.  The adsorption and competitive adsorption of bile salts and whey protein at the oil-water interface.  Soft Matter 2011, 7, 8942–8951.
  * Doudou et al. (2012) Doudou, S.; Vaughan, D.J.; Livens, F.R.; Burton, N.A.  Atomistic Simulations of Calcium Uranyl(VI) Carbonate Adsorption on Calcite and Stepped-Calcite Surfaces.  Environ. Sci. Tech. 2012, 46, 7587–7594.
  * Pomes and Roux (1998) Pomes, R.; Roux, B.  Free energy profiles for H+ conduction along hydrogen-bonded chains of water molecules.  Biophys. J. 1998, 75, 33–40.
  * Jagoda-Cwiklik et al. (2011) Jagoda-Cwiklik, B.; Cwiklik, L.; Jungwirth, P.  Behavior of the Eigen Form of Hydronium at the Air/Water Interface.  J. Phys. Chem. A 2011, 115, 5881–5886.
  * Calvo and Mottet (2011) Calvo, F.; Mottet, C.  Order-disorder transition in Co-Pt nanoparticles: Coexistence, transition states, and finite-size effects.  Phys. Rev. B 2011, 84.
  * Sharma and Debenedetti (2012) Sharma, S.; Debenedetti, P.G.  Free Energy Barriers to Evaporation of Water in Hydrophobic Confinement.  J. Phys. Chem. B 2012, 116, 13282–13289.
  * Ridder et al. (2002) Ridder, L.; Rietjens, I.; Vervoort, J.; Mulholland, A.  Quantum mechanical/molecular mechanical free energy Simulations of the glutathione S-transferase (M1-1) reaction with phenanthrene 9,10-oxide.  J. Am. Chem. Soc. 2002, 124, 9926–9936.
  * Kaestner et al. (2006) Kaestner, J.; Senn, H.; Thiel, S.; Otte, N.; Thiel, W.  QM/MM free-energy perturbation compared to thermodynamic integration and umbrella sampling: Application to an enzymatic reaction.  J. Chem. Theory Comput. 2006, 2, 452–461.
  * Wang et al. (2007) Wang, S.; Hu, P.; Zhang, Y.  Ab initio quantum mechanical/molecular mechanical molecular dynamics simulation of enzyme catalysis: The case of histone lysine methyltransferase SET7/9.  J. Phys. Chem. B 2007, 111, 3758–3764.
  * Ke et al. (2011) Ke, Z.; Guo, H.; Xie, D.; Wang, S.; Zhang, Y.  Ab Initio QM/MM Free-Energy Studies of Arginine Deiminase Catalysis: The Protonation State of the Cys Nucleophile.  J. Phys. Chem. B 2011, 115, 3725–3733.
  * Yan et al. (2011) Yan, S.; Li, T.; Yao, L.  Mutational Effects on the Catalytic Mechanism of Cellobiohydrolase I from Trichoderma reesei.  J. Phys. Chem. B 2011, 115, 4982–4989.
  * Mujika et al. (2012) Mujika, J.I.; Lopez, X.; Mulholland, A.J.  Mechanism of C-terminal intein cleavage in protein splicing from QM/MM molecular dynamics simulations.  Org. Biomolec. Chem. 2012, 10, 1207–1218.
  * Lonsdale et al. (2012) Lonsdale, R.; Hoyle, S.; Grey, D.T.; Ridder, L.; Mulholland, A.J.  Determinants of Reactivity and Selectivity in Soluble Epoxide Hydrolase from Quantum Mechanics/Molecular Mechanics Modeling.  Biochemistry 2012, 51, 1774–1786.
  * Rooklin et al. (2012) Rooklin, D.W.; Lu, M.; Zhang, Y.  Revelation of a Catalytic Calcium-Binding Site Elucidates Unusual Metal Dependence of a Human Apyrase.  J. Am. Chem. Soc. 2012, 134, 15595–15603.
  * Lior-Hoffmann et al. (2012) Lior-Hoffmann, L.; Wang, L.; Wang, S.; Geacintov, N.E.; Broyde, S.; Zhang, Y.  Preferred WMSA catalytic mechanism of the nucleotidyl transfer reactionin human DNA polymerase kappa elucidates error-free bypass of a bulky DNA lesion.  Nucl. Acids Res. 2012, 40, 9193–9205.
  * Crouzy et al. (2001) Crouzy, S.; Berneche, S.; Roux, B.  Extracellular blockade of K+ channels by TEA: Results from molecular dynamics simulations of the KcsA channel.  J. Gen. Physiol. 2001, 118, 207–217.
  * Allen et al. (2003) Allen, T.; Bastug, T.; Kuyucak, S.; Chung, S.  Gramicidin A channel as a test ground for molecular dynamics force fields.  Biophys. J. 2003, 84, 2159–2168.
  * Hub and De Groot (2008) Hub, J.S.; De Groot, B.L.  Mechanism of selectivity in aquaporins and aquaglyceroporins.  Proc. Natl. Acad. Sci. USA 2008, 105, 1198–1203.
  * Xin et al. (2011) Xin, L.; Su, H.; Nielsen, C.H.; Tang, C.; Torres, J.; Mu, Y.  Water permeation dynamics of AqpZ: A tale of two states.  BBA-Biomembranes 2011, 1808, 1581–1586.
  * Furini and Domene (2011) Furini, S.; Domene, C.  Selectivity and Permeation of Alkali Metal Ions in K+-channels.  J. Mol. Biol. 2011, 409, 867–878.
  * Kim and Allen (2011) Kim, I.; Allen, T.W.  On the selective ion binding hypothesis for potassium channels.  Proc. Natl. Acad. Sci. USA 2011, 108, 17963–17968.
  * Domene and Furini (2012) Domene, C.; Furini, S.  Molecular Dynamics Simulations of the TrkH Membrane Protein.  Biochemistry 2012, 51, 1559–1565.
  * Zhu and Hummer (2012) Zhu, F.; Hummer, G.  Theory and Simulation of Ion Conduction in the Pentameric GLIC Channel.  J. Chem. Theor. Comput. 2012, 8, 3759–3768.
  * Zhongjin and Jian (2011) Zhongjin, H.; Jian, Z.  Steered Molecular Dynamics Simulations of Ions Traversing Through Carbon Nanotubes.  Acta Chim. Sin. 2011, 69, 2901–2907.
  * Nalaparaju and Jiang (2012) Nalaparaju, A.; Jiang, J.  Ion Exchange in Metal-Organic Framework for Water Purification: Insight from Molecular Simulation.  J. Phys. Chem. C 2012, 116, 6925–6931.
  * Barducci et al. (2011) Barducci, A.; Bonomi, M.; Parrinello, M.  Metadynamics.  Wires. Comput. Mol. Sci. 2011, 1, 826–843.
  * Huber et al. (1994) Huber, T.; Torda, A.; van Gunsteren, W.F.  Local elevation: A method for improving the searching properties of molecular dynamics simulation.  J. Comput.-Aided Mol. Des. 1994, 8, 695–708.
  * Laio et al. (2005) Laio, A.; Rodriguez-Fortea, A.; Gervasio, F.L.; Ceccarelli, M.; Parrinello, M.  Assessing the accuracy of metadynamics.  J. Phys. Chem. B 2005, 109, 6714–6721.
  * Bussi et al. (2006) Bussi, G.; Laio, A.; Parrinello, M.  Equilibrium free energies from nonequilibrium metadynamics.  Phys. Rev. Lett. 2006, 96, 090601.
  * Raiteri et al. (2006) Raiteri, P.; Laio, A.; Gervasio, F.L.; Micheletti, C.; Parrinello, M.  Efficient reconstruction of complex free energy landscapes by multiple walkers metadynamics.  J. Phys. Chem. B 2006, 110, 3533–3539.
  * Barducci et al. (2008) Barducci, A.; Bussi, G.; Parrinello, M.  Well-tempered metadynamics: A smoothly converging and tunable free-energy method.  Phys. Rev. Lett. 2008, 100, 020603.
  * Singh et al. (2011) Singh, S.; Chiu, C.c.; de Pablo, J.J.  Flux Tempered Metadynamics.  J. Stat. Phys. 2011, 145, 932–945.
  * Bonomi et al. (2009) Bonomi, M.; Barducci, A.; Parrinello, M.  Reconstructing the Equilibrium Boltzmann Distribution from Well-Tempered Metadynamics.  J. Comput. Chem. 2009, 30, 1615–1621.
  * Branduardi et al. (2012) Branduardi, D.; Bussi, G.; Parrinello, M.  Metadynamics with adaptive Gaussians.  J. Chem. Theory Comput. 2012, 8, 2247–2254.
  * McGrath et al. (2013) McGrath, M.J.; Kuo, I.F.W.; Hayashi, S.; Takada, S.  ATP hydrolysis mechanism in kinesin studied by combined quantum-mechanical/molecular-mechanical metadynamics simulations.  J. Am. Chem. Soc. 2013.
  * Mantz et al. (2009) Mantz, Y.A.; Branduardi, D.; Bussi, G.; Parrinello, M.  Ensemble of Transition State Structures for the Cis- Trans Isomerization of N-Methylacetamide.  J. Phys. Chem. B 2009, 113, 12521–12529.
  * Leone et al. (2009) Leone, V.; Lattanzi, G.; Molteni, C.; Carloni, P.  Mechanism of action of cyclophilin a explored by metadynamics simulations.  PLOS Comput. Biol. 2009, 5, e1000309.
  * Melis et al. (2009) Melis, C.; Bussi, G.; Lummis, S.C.; Molteni, C.  Trans- cis Switching Mechanisms in Proline Analogues and Their Relevance for the Gating of the 5-HT3 Receptor.  J. Phys. Chem. B 2009, 113, 12148–12153.
  * Bussi et al. (2006) Bussi, G.; Gervasio, F.L.; Laio, A.; Parrinello, M.  Free-energy landscape for $\beta$ hairpin folding from combined parallel tempering and metadynamics.  J. Am. Chem. Soc. 2006, 128, 13435–13441.
  * Gangupomu and Abrams (2010) Gangupomu, V.; Abrams, C.  All-atom models of the membrane-spanning domain of HIV-1 gp41 from metadynamics.  Biophys. J. 2010, 99, 3438–3444.
  * Berteotti et al. (2011) Berteotti, A.; Barducci, A.; Parrinello, M.  Effect of Urea on the $\beta$-Hairpin Conformational Ensemble and Protein Denaturation Mechanism.  J. Am. Chem. Soc. 2011, 133, 17200–17206.
  * Granata et al. (2013) Granata, D.; Camilloni, C.; Vendruscolo, M.; Laio, A.  Characterization of the free-energy landscapes of proteins by NMR-guided metadynamics.  Proc. Natl. Acad. Sci. USA 2013, 110, 6817–6822.
  * Baftizadeh et al. (2012) Baftizadeh, F.; Biarnes, X.; Pietrucci, F.; Affinito, F.; Laio, A.  Multidimensional View of Amyloid Fibril Nucleation in Atomistic Detail.  J. Am. Chem. Soc. 2012, 134, 3886–3894.
  * Gervasio et al. (2005) Gervasio, F.L.; Laio, A.; Parrinello, M.  Flexible docking in solution using metadynamics.  J. Am. Chem. Soc. 2005, 127, 2600–2607.
  * Soederhjelm et al. (2012) Soederhjelm, P.; Tribello, G.A.; Parrinello, M.  Locating binding poses in protein-ligand systems using reconnaissance metadynamics.  Proc. Natl. Acad. Sci. USA 2012, 109, 5170–5175.
  * Limongelli et al. (2013) Limongelli, V.; Bonomi, M.; Parrinello, M.  Funnel metadynamics as accurate binding free-energy method.  Proc. Natl. Acad. Sci. USA 2013, 110, 6358–6363.
  * Sutto and Gervasio (2013 (published ahead of print) Sutto, L.; Gervasio, F.L.  Effects of oncogenic mutations on the conformational free-energy landscape of EGFR kinase.  Proc. Natl. Acad. Sci. USA 2013 (published ahead of print).  DOI: 10.1073/pnas.1221953110.
  * Martonak et al. (2006) Martonak, R.; Donadio, D.; Oganov, A.R.; Parrinello, M.  Crystal structure transformations in SiO2 from classical and ab initio metadynamics.  Nature Mat. 2006, 5, 623–626.
  * Trudu et al. (2006) Trudu, F.; Donadio, D.; Parrinello, M.  Freezing of a Lennard-Jones fluid: from nucleation to spinodal regime.  Phys. Rev. Lett. 2006, 97, 105701.
  * Stack et al. (2011) Stack, A.G.; Raiteri, P.; Gale, J.D.  Accurate rates of the complex mechanisms for growth and dissolution of minerals using a combination of rare-event theories.  J. Am. Chem. Soc. 2011, 134, 11–14.
  * Zhang et al. (2012) Zhang, C.; Knyazev, D.G.; Vereshaga, Y.A.; Ippoliti, E.; Nguyen, T.H.; Carloni, P.; Pohl, P.  Water at hydrophobic interfaces delays proton surface-to-bulk transfer and provides a pathway for lateral proton diffusion.  Proc. Natl. Acad. Sci. USA 2012, 109, 9744–9749.
  * Das et al. (2006) Das, P.; Moll, M.; Stamati, H.; Kavraki, L.E.; Clementi, C.  Low-dimensional, free-energy landscapes of protein-folding reactions by nonlinear dimensionality reduction.  Proc. Natl. Acad. Sci. USA 2006, 103, 9885–9890.
  * Perilla and Woolf (2012) Perilla, J.R.; Woolf, T.B.  Towards the prediction of order parameters from molecular dynamics simulations in proteins.  J. Chem. Phys. 2012, 136.
  * Ceriotti et al. (2011) Ceriotti, M.; Tribello, G.A.; Parrinello, M.  Simplifying the representation of complex free-energy landscapes using sketch-map.  Proc. Natl. Acad. Sci. USA 2011, 108, 13023–13028.
  * Tribello et al. (2010) Tribello, G.A.; Ceriotti, M.; Parrinello, M.  A self-learning algorithm for biased molecular dynamics.  Proc. Natl. Acad. Sci. USA 2010, 107, 17509–17514.
  * Tribello et al. (2011) Tribello, G.A.; Cuny, J.; Eshet, H.; Parrinello, M.  Exploring the free energy surfaces of clusters using reconnaissance metadynamics.  J. Chem. Phys. 2011, 135.
  * Bartels and Karplus (1998) Bartels, C.; Karplus, M.  Probability distributions for complex systems: adaptive umbrella sampling of the potential energy.  J. Phys. Chem. B 1998, 102, 865–880.
  * Micheletti et al. (2004) Micheletti, C.; Laio, A.; Parrinello, M.  Reconstructing the density of states by history-dependent metadynamics.  Phys. Rev. Lett. 2004, 92, 170601.
  * Bonomi and Parrinello (2010) Bonomi, M.; Parrinello, M.  Enhanced sampling in the well-tempered ensemble.  Phys. Rev. Lett. 2010, 104, 190601.
  * Do et al. (2013) Do, T.N.; Carloni, P.; Varani, G.; Bussi, G.  RNA/Peptide Binding Driven by Electrostatic Insight from Bidirectional Pulling Simulations.  J. Chem. Theory Comput. 2013, 9, 1720–1730.
  * Roitberg and Elber (1991) Roitberg, A.; Elber, R.  Modeling side-chains in peptides and proteins - application of the locally enhanced sampling and the simulated annealing methods to find minimum energy conformations.  J. Chem. Phys. 1991, 95, 9277–9287.
  * Patel et al. (2011) Patel, A.J.; Varilly, P.; Chandler, D.; Garde, S.  Quantifying Density Fluctuations in Volumes of All Shapes and Sizes Using Indirect Umbrella Sampling.  J. Stat. Phys. 2011, 145, 265–275.
  * Mueller et al. (2012) Mueller, M.; Smirnova, Y.G.; Marelli, G.; Fuhrmans, M.; Shi, A.C.  Transition Path from Two Apposed Membranes to a Stalk Obtained by a Combination of Particle Simulations and String Method.  Phys. Rev. Lett. 2012, 108.
  * Pietrucci and Laio (2009) Pietrucci, F.; Laio, A.  A collective variable for the efficient exploration of protein beta-sheet structures: application to sh3 and gb1.  J. Chem. Theory Comput. 2009, 5, 2197–2201.
  * Branduardi et al. (2007) Branduardi, D.; Gervasio, F.L.; Parrinello, M.  From A to B in free energy space.  J. Chem. Phys. 2007, 126, 054103.
  * Zinovjev et al. (2012) Zinovjev, K.; Martí, S.; Tuñón, I.  A Collective Coordinate to Obtain Free Energy Profiles for Complex Reactions in Condensed Phases.  J. Chem. Theory Comput. 2012, 8, 1795–1801.
  * Spiwok and Králová (2011) Spiwok, V.; Králová, B.  Metadynamics in the conformational space nonlinearly dimensionally reduced by Isomap.  J. Chem. Phys. 2011, 135, 224504–224504.
  * Kirkpatrick et al. (1983) Kirkpatrick, S.; Jr., D.G.; Vecchi, M.P.  Optimization by simulated annealing.  Science 1983, 220, 671–680.
  * Marinari and Parisi (1992) Marinari, E.; Parisi, G.  Simulated tempering: a new Monte Carlo scheme.  Europhys. Lett. 1992, 19, 451.
  * Park and Pande (2007) Park, S.; Pande, V.S.  Choosing weights for simulated tempering.  Phys. Rev. E 2007, 76, 016703.
  * Hansmann (1997) Hansmann, U.H.  Parallel tempering algorithm for conformational studies of biological molecules.  Chem. Phys. Lett. 1997, 281, 140–150.
  * Sugita and Okamoto (1999) Sugita, Y.; Okamoto, Y.  Replica-exchange molecular dynamics method for protein folding.  Chem. Phys. Lett. 1999, 314, 141–151.
  * Frenkel (2004) Frenkel, D.  Speed-up of Monte Carlo simulations by sampling of rejected states.  Proc. Natl. Acad. Sci. USA 2004, 101, 17571–17575.
  * Coluzza and Frenkel (2005) Coluzza, I.; Frenkel, D.  Virtual-Move Parallel Tempering.  ChemPhysChem 2005, 6, 1779–1783.
  * Sindhikara et al. (2008) Sindhikara, D.; Meng, Y.; Roitberg, A.E.  Exchange frequency in replica exchange molecular dynamics.  J. Chem. Phys. 2008, 128, 024103.
  * Bussi (2009) Bussi, G.  A simple asynchronous replica-exchange implementation.  Nuovo Cimento C Geophysics Space Physics C 2009, 32, 61–65.
  * Gallicchio et al. (2008) Gallicchio, E.; Levy, R.M.; Parashar, M.  Asynchronous replica exchange for molecular simulations.  J. Comput. Chem. 2008, 29, 788–794.
  * Rosta et al. (2009) Rosta, E.; Buchete, N.V.; Hummer, G.  Thermostat artifacts in replica exchange molecular dynamics simulations.  J. Chem. Theory Comput. 2009, 5, 1393–1399.
  * Sindhikara et al. (2010) Sindhikara, D.J.; Emerson, D.J.; Roitberg, A.E.  Exchange often and properly in replica exchange molecular dynamics.  J. Chem. Theory Comput. 2010, 6, 2804–2808.
  * Vreede et al. (2005) Vreede, J.; Crielaard, W.; Hellingwerf, K.; Bolhuis, P.  Predicting the signaling state of photoactive yellow protein.  Biophys. J. 2005, 88, 3525–3535.
  * Zhang and Mu (2012) Zhang, T.; Mu, Y.  Initial Binding of Ions to the Interhelical Loops of Divalent Ion Transporter CorA: Replica Exchange Molecular Dynamics Simulation Study.  PLOS ONE 2012, 7, e43872.
  * Zhou (2003) Zhou, R.  Trp-cage: Folding free energy landscape in explicit water.  Proc. Natl. Acad. Sci. USA 2003, 100, 13280–13285.
  * Garcia and Onuchic (2003) Garcia, A.; Onuchic, J.  Folding a protein in a computer: An atomic description of the folding/unfolding of protein A.  Proc. Natl. Acad. Sci. USA 2003, 100, 13898–13903.
  * Im et al. (2003) Im, W.; Feig, M.; Brooks, C.  An implicit membrane generalized born theory for the study of structure, stability, and interactions of membrane proteins.  Biophys. J. 2003, 85, 2900–2918.
  * Mei et al. (2012) Mei, Y.; Wei, C.; Yip, Y.M.; Ho, C.Y.; Zhang, J.Z.H.; Zhang, D.  Folding and thermodynamic studies of Trp-cage based on polarized force field.  Theo. Chem. Accts. 2012, 131.
  * Berhanu et al. (2013) Berhanu, W.M.; Jiang, P.; Hansmann, U.H.E.  Folding and association of a homotetrameric protein complex in an all-atom Go model.  Phys. Rev. E 2013, 87.
  * Kokubo and Okamoto (2004) Kokubo, H.; Okamoto, Y.  Self-assembly of transmembrane helices of bacteriorhodopsin by a replica-exchange Monte Carlo simulation.  Chem. Phys. Lett. 2004, 392, 168–175.
  * Oshaben et al. (2012) Oshaben, K.M.; Salari, R.; McCaslin, D.R.; Chong, L.T.; Horne, W.S.  The Native GCN4 Leucine-Zipper Domain Does Not Uniquely Specify a Dimeric Oligomerization State.  Biochemistry 2012, 51, 9581–9591.
  * Sugita and Okamoto (2000) 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, 261–270.
  * Fukunishi et al. (2002) Fukunishi, H.; Watanabe, O.; Takada, S.  On the Hamiltonian replica exchange method for efficient sampling of biomolecular systems: Application to protein structure prediction.  J. Chem. Phys. 2002, 116, 9058.
  * Liu et al. (2005) Liu, P.; Kim, B.; Friesner, R.A.; Berne, B.  Replica exchange with solute tempering: A method for sampling biological systems in explicit water.  Proc. Natl. Acad. Sci. USA 2005, 102, 13749–13754.
  * Affentranger et al. (2006) Affentranger, R.; Tavernelli, I.; Di Iorio, E.E.  A novel Hamiltonian replica exchange MD protocol to enhance protein conformational space sampling.  J. Chem. Theory Comput. 2006, 2, 217–228.
  * Fajer et al. (2008) Fajer, M.; Hamelberg, D.; McCammon, J.A.  Replica-exchange accelerated molecular dynamics (REXAMD) applied to thermodynamic integration.  J. Chem. Theory Comput. 2008, 4, 1565–1569.
  * Xu et al. (2008) Xu, C.; Wang, J.; Liu, H.  A Hamiltonian Replica Exchange Approach and Its Application to the Study of Side-Chain Type and Neighbor Effects on Peptide Backbone Conformations.  J. Chem. Theory Comput. 2008, 4, 1348–1359.
  * Zacharias (2008) Zacharias, M.  Combining elastic network analysis and molecular dynamics simulations by hamiltonian replica exchange.  J. Chem. Theory Comput. 2008, 4, 477–487.
  * Vreede et al. (2009) Vreede, J.; Wolf, M.G.; de Leeuw, S.W.; Bolhuis, P.G.  Reordering hydrogen bonds using Hamiltonian replica exchange enhances sampling of conformational changes in biomolecular systems.  J. Phys. Chem. B 2009, 113, 6484–6494.
  * Itoh et al. (2010) Itoh, S.G.; Okumura, H.; Okamoto, Y.  Replica-exchange method in van der Waals radius space: Overcoming steric restrictions for biomolecules.  J. Chem. Phys. 2010, 132, 134105.
  * Meng and Roitberg (2010) Meng, Y.; Roitberg, A.E.  Constant pH replica exchange molecular dynamics in biomolecules using a discrete protonation model.  J. Chem. Theory Comput. 2010, 6, 1401–1412.
  * Terakawa et al. (2011) Terakawa, T.; Kameda, T.; Takada, S.  On easy implementation of a variant of the replica exchange with solute tempering in GROMACS.  J. Comput. Chem. 2011, 32, 1228–1234.
  * Wang et al. (2011) Wang, L.; Friesner, R.A.; Berne, B.  Replica exchange with solute scaling: A more efficient version of replica exchange with solute tempering (REST2).  J. Phys. Chem. B 2011, 115, 9431–9438.
  * Zhang and Ma (2012) Zhang, C.; Ma, J.  Folding helical proteins in explicit solvent using dihedral-biased tempering.  Proc. Natl. Acad. Sci. USA 2012, 109, 8139–8144.
  * Bussi (2013) Bussi, G.  Hamiltonian replica-exchange in GROMACS: a flexible implementation.  Mol. Phys. 2013.  DOI:10.1080/00268976.2013.824126.
  * Huang et al. (2007) Huang, X.; Hagen, M.; Kim, B.; Friesner, R.A.; Zhou, R.; Berne, B.  Replica exchange with solute tempering: efficiency in large scale systems.  J. Phys. Chem. B 2007, 111, 5405–5410.
  * Denschlag et al. (2009) Denschlag, R.; Lingenheil, M.; Tavan, P.; Mathias, G.  Simulated solute tempering.  J. Chem. Theory Comput. 2009, 5, 2847–2857.
  * Zuckerman and Lyman (2006) Zuckerman, D.M.; Lyman, E.  A second look at canonical sampling of biomolecules using replica exchange simulation.  J. Chem. Theory Comput. 2006, 2, 1200–1202.
  * Chodera et al. (2007) Chodera, J.D.; Swope, W.C.; Pitera, J.W.; Seok, C.; Dill, K.A.  Use of the weighted histogram analysis method for the analysis of simulated and parallel tempering simulations.  J. Chem. Theory Comput. 2007, 3, 26–41.
  * Camilloni et al. (2008) Camilloni, C.; Provasi, D.; Tiana, G.; Broglia, R.A.  Exploring the protein G helix free-energy surface by solute tempering metadynamics.  Proteins 2008, 71, 1647–1654.
  * Deighan et al. (2012) Deighan, M.; Bonomi, M.; Pfaendtner, J.  Efficient Simulation of Explicitly Solvated Proteins in the Well-Tempered Ensemble.  J. Chem. Theory Comput. 2012, 8, 2189–2192.
  * Piana and Laio (2007) Piana, S.; Laio, A.  A bias-exchange approach to protein folding.  J. Phys. Chem. B 2007, 111, 4553–4559.
  * Baftizadeh et al. (2012) Baftizadeh, F.; Cossio, P.; Pietrucci, F.; Laio, A.  Protein Folding and Ligand-Enzyme Binding from Bias-Exchange Metadynamics Simulations.  Curr. Phys. Chem. 2012, 2, 79–91.
  * E et al. (2002) E, W.; Ren, W.; Vanden-Eijnden, E.  String method for the study of rare events.  Phys. Rev. B. 2002, 66, 052301.
  * Maragliano et al. (2006) Maragliano, L.; Fischer, A.; Vanden-Eijnden, E.; Ciccotti, G.  String method in collective variables: Minimum free energy paths and isocommittor surfaces.  J. Chem. Phys. 2006, 125, 024106.
  * Vashisth and Abrams (2013) Vashisth, H.; Abrams, C.F.  All-atom structural models of insulin binding to the insulin receptor in the presence of a tandem hormone-binding element.  Proteins 2013, 81, 1017–1030.