although 
 in principle observing nature carefully over infinite 
 time scales could be sufficient to reveal all physics , it is more 
 practical to design experiments that investigate specific questions 
 by intentionally varying physical parameters in a controlled manner . 

similarly in computational modeling , when direct simulation of natural 
 processes by molecular dynamics is infeasible , specially designed 
 simulations can nonetheless reveal key physics in far less simulation 
 time .
the adaptive enhanced sampling method metadynamics is one such approach specifically designed for the determination 
 of potentials of mean force ( pmfs ) by promoting transitions between 
 long - lived metastable states .
metadynamics has been widely applied 
 across chemistry from materials science to biochemistry , but it is 
 still quite young theoretically , with a rigorous proof of convergence 
 appearing only one year ago .
metadynamics 
 works by using a choice of reduced coordinates called 
 collective variables ( cvs ) to iteratively build a bias that increases 
 the rates of transitions between metastable energy wells ; increased 
 transition rates imply decreased sampling autocorrelation and thus 
 improved pmf estimates .
other adaptive methods of the same generation 
 and similar philosophy include the adaptive biasing force and wang  landau algorithms , 
 and it is descended from the local elevation method . the degree to which a bias can actually promote those transitions , 
 however , depends on how well the cvs capture the true reaction coordinates . 

when the cvs are imperfect , the results may not approach the true 
 pmf rapidly , and this is common enough that it is often considered 
 to be the single most relevant limitation preventing application of 
 metadynamics to the study of complex systems .
furthermore , the cost of building the bias also depends on the complexity 
 of the cvs , scaling with the volume of cv space 
all enhanced sampling methods that rely on cvs share 
 these drawbacks to greater or lesser extents .
this paper describes a new variant of metadynamics , 
 metabasin metadynamics ( mbmetad ) , that is designed to suffer less 
 from the use of poor quality cvs and to make more practical the use 
 of larger numbers of cvs by judiciously restricting the bias s 
 domain in cv space .
however , in order to discuss the features of metadynamics 
 that we wish to improve with our new method , we must first compare 
 to a more venerable alternative , window - based umbrella sampling .
it is accomplished by running 
 many simulations with different energetic biases that keep each simulation 
 restrained within a different small region , or window , of cv space . 

these windows are constructed so that the sampled distributions in 
 the windows overlap with one another to cover all of the phase space 
 of interest in a given investigation ; this typically involves choosing 
 1 ) a scale for each stratified dimension of cv space to set the separation 
 of window centers from one another and 2 ) a single energy scale to 
 set the strengths of the restraint .
once these are chosen , one then 
 runs simulations in each window , perhaps with some form of replica 
 exchange among the biased walkers .
after this , 
 the sampling in all of the windows forms a patchwork covering of the 
 cv space that can then be sewn together into a single overall pmf 
 estimate using a method such as gaussian process regression , the weighted histogram analysis method , or the multistate bennet acceptance ratio .
metadynamics , on the other hand , is an 
 auxiliary distribution sampling 
 approach .
it functions by iteratively building a bias away from previously 
 visited points to accelerate escape from metastable basins and thereby 
 decrease autocorrelation of sampling on the cv space . in tempered 
 metadynamics ,
the type of metadynamics that 
 converges asymptotically , the bias is 
 constructed by adding gaussian hills of bias of progressively shrinking 
 height at regular intervals .
overall , a typical case requires choosing 
 1 ) the length scale of the hill in each dimension , 2 ) a single hill 
 energy scale , and 3 ) a single energy scale for the convergence rate 
 of the bias .
after these are chosen , one then runs one or more simulations 
 under an iteratively growing bias , again with the possibility of replica 
 exchange .
once this is done , one can calculate pmfs via their 
 direct connection to the bias in the asymptotic regime or via one 
 of several nonequilibrium reweighting estimators .
when 
 stratifying based on a set of cvs is not enough to ensure 
 good sampling over the cv range , one says that there are hidden slow 
 variables ( hsvs ) .
when these are present , replica exchange can be 
 inefficient , and sampling in the windows often does not accurately 
 represent true boltzmann sampling on accessible simulation time scales . the free energy differences between windows spaced far apart are 
 pieced together using information from windows between them , making 
 the errors induced by poor sampling in a single window nonlocal .
this 
 nonlocality follows the topology of the pmf in a way such that sampling 
 errors in transition regions have large effects on relative free energy 
 estimates between the basins they connect , while errors in 
 basin regions and high energy regions have a more local effect .
however , unless the hsvs prevent replica 
 exchange or lead to obvious unphysical features in a pmf , the erroneous 
 sampling can easily be mistaken for correct .
the presence of 
 hsvs in metadynamics , on the other hand , typically 
 causes sampling hysteresis that is clearly visible in the bias .
when the cvs are open for modification , 
 this is an advantage because it gives a clear signal that the cvs 
 need improvement . as in window - based umbrella sampling ,
the effect 
 on pmf estimates is nonlocal , but , unlike in that method , the nonlocal effects from hsv sampling in high - energy regions can 
 affect the estimates of the pmf across all other regions . 

specifically , the nonequilibrium bias can drive undesired changes 
 in hsvs that never relax on the time scale of simulation , for instance 
 driving undesired boiling of solvent or an irreversible fluctuation 
 of a protein conformation .
therefore , although hysteresis in the presence of 
 hsvs is a desired behavior in regions of cv space relevant to the 
 process of interest , it is an inconvenience in high - energy regions . 

it is especially troublesome in complex systems with sensitive dynamical 
 environments and many potential hsvs , e.g. , almost all proteins .
regardless of the presence of hsvs , the simplest form of windowed 
 umbrella sampling scales poorly with the dimensionality of cv space 
 because one typically covers the space with the fixed - volume windows , 
 and the cost of covering space with fixed - volume sets scales exponentially 
 with dimension .
this can be ameliorated somewhat by taking a metadynamics - like 
 approach in which one runs simulations only in windows in low - energy 
 regions encompassing the main basins of interest , i.e. , a metabasin , 
 gradually learning the shape of the metabasin along the way .
the metabasin manifold has a lower volume  and 
 often a lower effective dimensionality  than its bounding box , 
 so covering it with windows is less prohibitive than covering the 
 full span of the cv space .
no advance knowledge of the metabasin s 
 shape is required , only an energy scale of interest . in this 
 paper
however , our goal and approach differ . our primary goal with mbmetad 
 is to prevent driving in high - energy regions to prevent the excitation 
 of hsvs that irreversibly alter sampling .
greater efficiency with 
 many cvs is another benefit and an important one , but it is currently 
 a less pressing issue in the application of metadynamics than the 
 problem of irreversible driving .
the strategy differs from adaptive 
 umbrella sampling because metadynamics already acts to fill metabasins 
 from the lowest point out  an initialization logic  and 
 our change is therefore rather to cause metadynamics to drive escape 
 no
the most common approach to this termination 
 in the literature 
 has been to self - limit by slowing down the addition of bias , for instance 
 in self - healing umbrella sampling , flat 
 histogram metadynamics , and well - tempered 
 metadynamics .
however , though these can 
 be effective , in such methods the time scales of updating the bias 
 and the energy scale of the self - limiting are tightly coupled , which 
 can make it impossible to simultaneously choose acceptable values 
 for both in complex systems .
most often , one must simply choose to 
 add energy very slowly or , similarly , decrease the rate of energy 
 addition very rapidly to prevent spoiling the dynamics .
there 
 is therefore a need for a scheme that takes a different approach , 
 sampling a specific energy scale regardless of how quickly bias is 
 added , to empower researchers to make more practically efficient investigations 
 of sensitive systems with metadynamics . in this paper
, we solve 
 the longstanding problem of designing metadynamics 
 to flatten only a region of arbitrary shape , without restricting sampling 
 to that region , in order to design a metadynamics that flattens only 
 up to a chosen energy level and no further regardless of how quickly 
 new bias is added .
the approach we take is compatible with all tempering 
 strategies , bias exchange , multiple walkers , and more , functioning as a modular enhancement that can and should be used 
 together with other advances in metadynamics methodology .
the task 
 of filling a finite domain has been a pernicious challenge throughout 
 the history of metadynamics . when used without proper understanding 
 it can suffer from boundary artifacts .
solutions are available for 
 1d intervals , but the first simple , implemented 
 solutions for rectangular boxes with hard wall boundary conditions 
 in any number of dimensions were first published by mcgovern and de 
 pablo only two years ago .
we independently 
 derived the same corrections in the course of studying the convergence 
 of metadynamics , and here we extend the mcgovern 
de pablo approach 
 to time - dependent domains of arbitrary shape with open boundaries . 

our method consists of this generalized boundary - corrected metadynamics 
 together with adaptive rules for defining time - dependent domains that 
 represent the metabasins of genuine physical interest in practical 
 metadynamics applications .
these rules do not require any a priori estimate of barrier heights in the system .
section 2 details this new 
 method , section 3 presents 
 results and discussion of three example applications to biomolecular 
 simulation , and section 4 concludes the paper .
tempered 
 metadynamics is proven to converge like a specific quasiequilibrium 
 differential equation that depends only on the method and not on underlying 
 system dynamics . therefore , to describe 
 the design of the new approach that converges on a finite , open - boundary 
 domain
, we first define a new shape of hill that gives rise to a differential 
 equation that flattens the energy on any domain .
second , because the 
 domain of interest will only rarely be known ahead of time , we next 
 discuss how to adaptively define these domains over the course of 
 simulation using either a known target basin energy level or an adaptively 
 learned energy level corresponding to the minimal metabasin containing 
 several known points of interest . in the original metadynamics ,
the adaptive 
 bias is constructed as a sum of gaussian hills centered at previously 
 visited points in cv space , i.e. , in 1d1where h is an energy - valued 
 hill height ,  is a cv - valued hill width , and the sti are the 
 sequence of points in cv space visited at each ti . to ensure convergence one must sequentially 
 decrease the hill height , for instance using the well - tempered metadynamics 
 ( wtmetad ) rule2where t is a parameter that controls the rate of height decrease . a general 
 theoretical form for metadynamics that can be proven to converge is3where g(s , s )
is called the hill kernel and encapsulates 
 the hill shape , and w[v(s , tn ) ] is 
 called the global tempering rule and serves as a complement or substitute 
 to the wtmetad rule for decreasing the hill sizes as time goes on ; 
 see for example transition - tempered metadynamics ( ttmetad ) or the global tempering used in experiment directed 
 metadynamics .
the new method of this 
 paper is modular with respect to the various tempering rules and only 
 involves redefinition of g(s , s ) .
an illustration of a metabasin metadynamics hill function 
 ( top , eq 4 ) on the unit 
 circle near 
 the domain boundary ( black lines denote r = 1 ) and 
 its decomposition into a rescaled hill ( left , eq 7 ) and a low plateau outside of the domain 
 ( right , eq 9 ) .
the hill 
 is based on a gaussian of width 1/10 placed at r = 0.9 .
the effects of adding rescaled hill and exterior plateau 
 contributions evenly over all points on the unit circle are exact 
 complements , leading to bias updates that are flat everywhere when 
 sampling is flat on the unit circle . specifically , consider the case that one knows of some hill 
 kernel g0(s , s ) 
 that could flatten an entire cv space but would like to flatten a 
 domain d of that cv space instead .
then the mbmetad 
 rule is to use a hill with two parts4where gint(s , s ) are hills 
 that would flatten the interior of the domain if it had hard walls , 
 and gext(s , s ) raises the bias level of the exterior 
 in such a way to exactly match the exterior bias level to the bias 
 level of the domain boundary while also exactly counteracting the 
 parts of gint(s , s ) that slope into the exterior 
 of the domain . both are based only on g0(s , s ) and the domain d. these hills are illustrated in figure 1 .
for the first part , we use a close 
 relative to the multiplicative 
 mcgovern  de pablo rule5where i(s ; d , g0 ) is a boundary - normalized 
 integral of g0(s , s ) over d for its s argument , i.e.6with d the boundary 
 of domain d. this integral does not have a closed 
 form solution for domains of general shape , but for gaussian hills 
 on a rectangular domain it is given by the sum of error functions 
 provided in the paper of mcgovern and de pablo .
the rule eq 5 has the undesirable property that as s 
 approaches the boundary of d from the interior , the 
 values of g(s , s ) in the exterior diverge to infinity , making it unsuitable 
 for metadynamics on an open boundary domain .
this occurs because i(s ; d , g0 ) goes to zero quickly as s goes beyond 
 the boundary , so an expedient solution is to modify i(s ; d , g0 ) such that it falls to a fixed value greater than zero .
in mbmetad , 
 therefore , our form for the first term in eq 4 is7with i(s ; d , g0 ) as above and 
 in the implementation described in this paper f(x ) is the simple piecewise polynomial f(x ) = x for x  
 1 and f(x ) = x + 
 0.5(1x )   for x < 1 .
this choice of f is continuous and differentiable 
 even at x = 1 , ensuring that all bias forces will 
 be continuous .
i(s ; d , g0 ) is constant in the interior of 
 the domain and decays to zero outside of the domain , so dividing by f(i(s ; d , g0 ) ) has the effect of leaving hills 
 in the interior unchanged and causing them to bulge out of the domain 
 subtly more than they otherwise would when close to its boundary . using only that first part of the hill , the differential equation 
 from the convergence paper8would predict that the interior 
 of the region 
 would be flattened , but the bias inside would grow indefinitely , forming 
 a steep - walled , flat - topped mesa in the pmf .
this would drive all 
 sampling out of the metabasin of interest and is clearly undesirable . 

specifically , we use9which plateaus 
 in the exterior of the region 
 with a plateau height equal to the boundary average value of the first 
 part of the new bias .
because 
 its height is set to match the bias added to the boundary , it ensures 
 that the exterior bias level always matches the average bias level 
 of the domain boundary . in other words ,
this part of the bias update 
 pulls sampling into the domain everywhere along the boundary exactly 
 enough to cancel out how the first part of the hill would push sampling 
 out of the domain close to the new hill .
moreover , the plateau is 
 perfectly shaped so that flat quasiequilibrium sampling on the interior 
 of the domain leads to an exactly flat increase of the bias everywhere 
 in cv space .
finally , this term , eq 9 , is only nonzero when a hill is added close to the 
 boundary ( s is close to the boundary ) , so 
 mbmetad functions just like normal metadynamics away from those boundaries . while this scheme might appear technical at first glance , it has 
 a simple multiscale physical interpretation . on a coarse scale ,
we 
 look at the system through a two - state lens : either the sample is 
 inside or outside the domain . on a fine - grained scale ,
we look at 
 the points in cv space to a resolution comparable to the length scale 
  of the original hill .
the metabasin rule , then , corresponds 
 to looking at the cv space using a topology in which all the points 
 inside the domain are resolved to a fine - grained level , while all 
 of the points outside are collapsed into one generic state .
the first 
 part of the hill is added to a fine - grained neighborhood of the interior , 
 while the second is added to the entire coarse - grained exterior state  also 
 a single neighborhood , but using the multiscale topology described 
 above instead of the usual metric topology of the cv space .
it is often argued that bias hills should match the green s 
 functions for the dynamics , in at least a loose sense , and this is no exception .
the simplest way to understand this new 
 hill shape in those terms is to recognize that these hills are also 
 like green s functions , but where boundary conditions have 
 been adjusted so that the walls are only partially reflective and 
 any particles that exit the boundary from the interior must be balanced 
 by particles crossing the boundary from the exterior . in the exterior , 

the constant profile in the far field corresponds to a well - mixed , 
 fixed concentration boundary condition at infinity .
the second part 
 is thus related to the quasistationary distribution of a well - mixed 
 domain exterior and therefore retains a natural connection to dynamics . up to now
we have spoken in terms of a single connected domain , 
 but the development above applies to any disconnected domain just 
 as easily .
simply divide the domain d into its n components di and whenever a sample would be added inside a domain di , use the above definitions with di in the place of d. in the discussion of its physical meaning , consider a 2-state coarse - graining rather than a two - state coarse - graining ; 
 the new set of states is the cartesian product of all the sets of 
 two states ( interior and exterior ) per domain .
when one s 
 goal is to flatten a pmf in a domain without 
 causing discontinuities , the only domains that can be flattened completely 
 are those with an isoenergetic boundary .
therefore , it is natural 
 to define valid domains by defining a boundary free energy and letting 
 the domain be the set of all cv points with a lower free energy than 
 that boundary energy . since the zero of free energy is arbitrary , 
 these definitions must be made in terms of an energy difference rather 
 than an energy alone .
for instance , in this paper , we will show the 
 use of domains like  all points less than 45 kj above the minimum 
 free energy  , a known basin free energy level , and 
all 
 points less than 15 kj above the transition barrier between points 
 a and b  , an implicitly defined metabasin energy level . of 
 course
it is often 
 one of the things one would like to discover using metadynamics . for 
 that reason
, we also add an on - the - fly mechanism to discover di defined in terms of simpler 
 criteria like the examples above .
we solve this problem much 
 as in ttmetad , where we used the bias 
 as an on - the - fly free energy estimator .
however , 
 in this case we seek metabasins rather than barriers , and now the 
 bias will not properly approximate the free energy anywhere outside 
 of a previously defined domain .
therefore , we use the nonequilibrium 
 umbrella sampling estimator of branduardi , bussi , and parrinello to approximate the domain locations on - the - fly 
 rather than the bias alone .
this allows for the discovery of extensions 
 of a domain outside a previously chosen domain , possibly including 
 new disconnected components , regardless of initial conditions .
the 
 free energy estimates provide minima , transition barriers , and on - the - fly 
 region selections as described in the next paragraph .
the recipes 
 above have straightforward numerical implementations . i(s ; d , g0 ) can be calculated by numerical integration on any domain d simply by summing many hills g0(s , s ) evenly over d , and once it is available f(i(s ; d , g0 ) ) and i(s ; d , g0)/f(i(s ; d , g0 ) ) 
 are simple to calculate . the nonequilibrium umbrella sampling estimator 
 is evaluated using a running histogram accumulated during simulation 
 alongside the bias , free energy minima are trivial to calculate , and 
 the transition barriers are found via breadth - first path search as 
 in ttmetad .
these minima and barriers , 
 plus an offset , define the target region energy .
the region of interest 
 is then the total set of grid points under that free energy in the 
 current on - the - fly estimate .
that set of grid points is then split 
 into connected domains with a standard flood - fill connected components 
 labeling algorithm .
once these domains 
 are defined , we compute i(s ; d , g0 ) for each domain and 
 then continue the simulation .
our implementation updates the domain 
 at fixed time intervals with an option not to update the domain each 
 time but to instead also wait until the free energy estimate in the 
 exterior of the region has changed by some set amount . with tempering , 

the latter corresponds to a logarithmic update schedule that cuts 
 down on unnecessary computational overhead at late times when the 
 domain definition is stable .
each change in domain definition only 
 affects future hill addition , so the previously deposited bias continues 
 to be used unchanged .
it provides mbmetad 
 with a minimum based domain 
 or a transition barrier based domain ( with transition barriers defined 
 as described in the paper introducing ttmetad ) , using only one additional required input number , the energy offset 
 from the transition or minimum .
two further parameters to control 
 the update schedule and one controlling accuracy of the numerical 
 integration are also available for fine - tuning , but these have general 
 defaults and can be considered optional .
our purpose with mbmetad is to fundamentally 
 change the trade - offs 
 of metadynamics with respect to other methods by giving it a self - limiting 
 mechanism , and thus the remainder of this paper focuses on investigating 
 the self - limiting behavior in three examples .
first , we show that 
 it functions as intended all the way through pmf convergence for a 
 simple model biomolecule with imperfect cvs , alanine dipeptide , without 
 negatively affecting convergence efficiency or accuracy .
second , we 
 show that it correctly self - limits in a more complex example , monomeric 
 actin , where the self - limiting is used to achieve more stable and 
 reproducible exploration of conformational space . finally , we apply 
 it to a membrane transport protein , clc - ec1 , using mediocre cvs to 
 examine its failure modes and
in particular show that it retains metadynamics 
 attractive property of displaying hysteresis in the domain of interest 
 ( rather than harder to recognize errors ) even as it prevents the nonequilibrium 
 driving from spoiling the dynamics in other domains .
our focus 
 in these examples is the self - limiting behavior and reproducibility 
 rather than convergence efficiency or accuracy ; self - limiting behavior 
 aimed at improving reproducibility is the sole novel feature introduced 
 in mbmetad , and the method otherwise behaves much like whatever other 
 metadynamics method it augments . though the accuracy of metadynamics 
 in practice does deserve further study , it is not a focus of this 
 particular paper .
a standard first 
 test case for new metadynamics methods is blocked alanine dipeptide in vacuo with  and  dihedral cvs because the 
 long literature on metadynamics provides excellent guidance on transferring 
 knowledge gained with this model to understanding of the method in 
 more complex systems .
alanine dipeptide provides a simple case of a system with imperfect 
 separation of time scales between cv dynamics and other relaxations 
 and is thus a minimal test for any adaptive enhanced sampling method 
 that is meant to be robust to memory in the cv dynamics .
these 
 tests modeled alanine dipeptide using the amber03 force field and at a constant temperature of 300 k using 
 a langevin thermostat with a drag of 5 amu / ps .
the dynamics were integrated 
 in gromacs 4.6.1 patched with a customized version of plumed2 and using a stochastic dynamics leapfrog algorithm with 
 a time step of 2 fs , particle - mesh ewald summation , and shake constraints on all bonds .
each simulation began from the same c7eq conformation , but each set of 32 
 runs began from 32 distinct initial velocities pseudorandomly generated 
 from distinct random seeds and used 32 distinct langevin thermostat 
 seeds .
we compared wtmetad , ttmetad , and transition - tempered 
 mbmetad using 
 hill height , width , and rate parameters and tempering parameters as 
 in our previous work introducing ttmetad .
the height was 1.2 kj / mol , the width 0.35 radians , hills were deposited 
 every 120 fs , the ttmetad parameter t was 2kbt , 
 and the transition wells for ttmetad were ( 1.25 , 1.25 ) and 
 ( 1.0 , 1.25 ) in (  ,  ) cv space ; for wtmetad t was 4 kbt . when using metabasin metadynamics we compared 
 the minimum - based domain definition with relative energy levels of 
 35.0 and 42.5 kj / mol and transition - based domain definition with relative 
 energy levels of 7.5 and 15.0 kj / mol .
we chose these values in order 
 to roughly match the two methods based on the approximation that the 
 apparent transition barrier relative to the minimum in these cvs is 
 27.5 kj / mol ; this makes for clear and direct comparison between the 
 two different styles of domain specification .
the lower energy 
 values were chosen to illustrate how the two domain definitions differ 
 in behavior as the energy excess above the transition barrier becomes 
 low and the pmf is barely flattened over the barrier between the transition 
 wells .
bias - based free energy estimates averaged over 32 runs of 8 ns 
 for four different metadynamics choices designed to give different 
 levels of resolution of the high - energy regions . clockwise from top 

left : wtmetad , ttmetad , 15 kj / mol transition - referenced mbmetad , 42.5 
 kj / mol minimum - referenced mbmetad .
the estimates are 
 all but identical in low - energy regions despite stark differences 
 in the high - energy regions .
wtmetad is least controllable and gives 
 the most resolution of high energy regions , ttmetad is more controllable 
 but would give more and more resolution of high energy regions given 
 additional convergence time , and both forms of mbmetad are fully controllable , 
 converging without resolving the high energy regions .
figures 
 on the left show convergence results for the estimated free energy 
 difference between points ( 1.25 , 1.0 ) and ( 1.25 , 1.0 ) , 
 the two primary basins , while figures on the right show convergence 
 results for the estimated free energy difference between points ( 1.25 , 
 1.0 ) and ( 2.1 , 1.5 ) , the deeper basin and the lowest free 
 energy barrier point between the two basins .
figures on top show the 
 averaged free energy difference estimates as a function of time , while 
 the lower figures show the standard deviation of the estimates as 
 a function of time .
all methods are as described in the text , with 
 the statistics calculated over 32 runs per method .
the mbmetad methods 
 in these figures used energy offsets of 42.5 kj / mol ( minimum referenced ) 
 and 15 kj / mol ( transition referenced ) . for these free energy difference 
 estimates inside the metabasin domain , mbmetad incurs no convergence 
 accuracy penalty .
the results for higher 
 energy limits confirm that mbmetad with 
 either domain definition retains the accuracy and efficiency of ttmetad 
 in the low - energy regions it is designed to study while giving up 
 all resolution of the high - energy regions it is designed not to explore .
figure 2 shows the replicate - averaged 
 estimated pmfs of regular ttmetad and the two mbmetad variants after 
 8 ns of simulation for comparison ; it is evident that the low - lying 
 contours are essentially indistinguishable among the three methods , 
 while the high - energy contours present in the first are simply missing 
 in the latter two . also as expected , the contours outside of the domain 
 but close enough to be biased are somewhat random and are not optimized 
 over the course of convergence .
figure 3 compares the convergence rates for these methods on 
 the domain of interest in terms of trueness and precision .
the mbmetad methods in these figures used energy 
 offsets of 42.5 kj / mol ( high energy minimum referenced ) , 35 kj / mol 
 ( low energy minimum referenced ) , 15 kj / mol ( high energy transition 
 referenced ) , and 7.5 kj / mol ( low energy transition referenced ) .
the 
 trueness of the mbmetad methods are all fairly consistent , but the 
 precision of transition - referenced mbmetad is more robust with respect 
 to choice of domain specification due to the early time cancellation 
 of errors discussed in the text . however ,
with lower energy limits , the transition - referenced 
 domain 
 definitions become clearly superior to minimum - referenced domains . 

this is because in the normal case that cvs are imperfect , hidden 
 barriers in hsvs can cause the bias to grow too high in one basin 
 before exploring the next basin .
when 
 this occurs in minimum - referenced mbmetad , the energy differences 
 from the minimum can therefore appear higher than they should and 
 the estimated domains then become too small .
if the intended domain 
 was not much larger than the true metabasin of interest , the bias 
 may not fully flatten all the way between the states of interest until 
 the overestimation corrects itself . on the other hand , when similar 
 overestimation occurs in transition - referenced mbmetad , the transition 
 barrier estimate also overshoots the true value , so that the domain 
 is made larger by the same amount that the previous error would make 
 it smaller .
this systematic cancellation of errors ensures greater 
 stability and better flattening before convergence and , thus better , 
 less autocorrelated sampling and faster convergence .
we always recommend 
 using the transition - based rule for adaptive domain restriction when 
 seeking to connect basins . in these tests , the desired regions 
 of interest are single large 
 connected domains including both wells . in the high energy tests , 

the on - the - fly regions begin with large connected domains including 
 both wells and small disconnected satellite domains due to sampling 
 noise in the free energy estimate . as the simulation proceeds and 

the free energy estimate is refined , new satellites spring up , some 
 disappear , and some connect with and absorb into the main domain . 
 the behavior is the same in the low energy transition - referenced simulation 
 as well . however , in the low energy minimum - referenced simulations , 
 the regions often begin with one primary domain around the basin the 
 simulation is initialized in . in these runs ,
when a transition to 
 a new well occurs there is a delay before the umbrella sampling estimator 
 registers the discovery of a new low - energy region , but then a domain 
 begins to grow around the new basin . depending on vagaries of sampling , 
 these domains can either remain separate or join together before the 
 simulation terminates .
this reconnection is not efficient , so once 
 again , we always recommend using the transition - based rule for adaptive 
 domain restriction when seeking to connect basins .
actin is a key component 
 of the cytoskeleton , but the structural 
 bases and mechanisms of its regulation are only partially understood , 
 even at the monomer level .
the monomer is also called globular actin 
 or g - actin , and aspects of it have been studied using metadynamics and umbrella sampling .
however , g - actin 
 is a dynamically sensitive protein with a high degree of allostery 
 that can be difficult to design good cvs to investigate and can resist 
 enhanced sampling analysis .
schematic illustration 
 of the cvs used to investigate nucleotide - dependent 
 actin dynamics , consisting of the distance between subdomains 2 and 
 4 and the torsion angle between subdomains 2 , 1 , 3 , and 4 .
a previous paper used 
 umbrella sampling 
 to investigate nucleotide - dependent g - actin conformational distributions 
 measured in terms of an interdomain distance and torsion ( figure 5 ) but was unable 
 to reach convergence , presumably due to the presence of hsvs .
attempting 
 to use metadynamics with the same torsion and distance cvs showed 
 hysteresis , confirming the presence of hsvs , but this hysteresis typically 
 involved unexpected and apparently irreversible behavior . hypothesizing 
 that the cvs were adequate in low - energy regions but inadequate in 
 higher energy configurations
, we applied mbmetad to the problem in 
 order to sample the lowest - lying metastable basins well and learn 
 enough to refine the cvs before investigating further . to do 
 this
, we performed 8 simulations of wtmetad with and without 
 adaptive domain restriction for 50 ns each .
the simulations were set 
 up as in previous work , with proteins 
 simulated in namd patched with a customized 
 version of plumed2 using the charmm22/27 
 force field with cmap .
our cvs were similar 
 to those in the previous work , consisting of the 2 - 1 - 3 - 4 torsion between 
 the cores of the four actin subdomains  the flatness of the 
 monomer  and the distance between subdomains 2 and 4the 
 width of the nucleotide binding cleft .
however , for the sake of computational 
 efficiency the positions of these subdomains were calculated as the 
 centers of mass of only the backbone and c carbons 
 in the subdomain cores rather than the centers of mass of all of the 
 subdomain core atoms .
all simulations were initialized from a flattened 
 conformation like that of filamentous actin ( pdb i d 2zwh ) .
metadynamics added hills of initial height 0.004184 kj / mol 
 and widths 0.3 radians and 0.02 nm every 200 fs , with heights sequentially 
 adjusted according to the wtmetad rule with a bias factor of 10 ; the 
 mbmetad domain was defined as every point less than 25.104 kj / mol 
 above the current estimate of the free energy minimum and was updated 
 every 200 ps if the exterior bias had increased by 2.092 kj / mol since 
 the last update .
a minimum reference is appropriate because the simulations 
 do not study a particular transition but rather the fluctuations about 
 a single basin .
representative sampling from unrestricted wtmetad ( top 
 row ) and 
 well tempered mbmetad ( lower two rows ) with the same initial conditions 
 ( columns ) as a function of cleft width ( w ) and twist angle (  ) . 
 in the unrestricted case only bias pmf estimates
are shown because 
 sampling and bias are directly related ; for mbmetad both bias pmf 
 estimates ( middle row ) and sampling histograms ( bottom row ) are shown .
representative biases and sampling 
 histograms at the end of 50 
 ns are shown in figure 6 .
it is evident that the mbmetad rule successfully limits the bias 
 to a smaller metastable basin of cv space than the unrestricted wtmetad 
 does without artificially restricting the sampling to the same basin . 

they consist of single primary domains around the starting states 
 together with smaller disconnected domains that spring up wherever 
 the simulation occupies another state long enough . however , the primary 
 purpose of this test is to demonstrate greater simulation repeatability 
 for mbmetad compared to the wtmetad reference , investigated in figure 7 .
figure showing the disagreement 
 between the free energy estimates 
 of different runs of like type .
disagreement ( y axis ) 
 is measured as one minus the ratio of the area of the intersection 
 of the two regions predicted by two different runs to be below a given 
 free energy threshold ( x axis ) to the area of their 
 union .
dashed lines represent the disagreement as a function of free 
 energy threshold for individual pairs , while solid lines indicate 
 their averages .
the unrestricted wtmetad estimates disagree with one 
 another more than the mbmetad estimates at all free energy levels , 
 showing much improved repeatability for the latter .
this figure demonstrates that overlaps of the biases 
 of different 
 randomly initialized runs are greater in mbmetad than wtmetad , i.e. , 
 that the sampling is more similar across runs of mbmetad than wtmetad , 
 just as desired .
the disagreement between predictions of which regions 
 of cv space are below a threshold energy shown in this figure correspond 
 to one minus the ratio of the area of the intersection of the regions 
 to the area of their union ; expressed as an integral , this is 1  
  min(i(fi(s ) < fth ) , i(fj(s ) < fth))ds/max(i(fi(s ) <
fth ) , i(fj(s ) < fth))ds , with fi and fj being the pmf estimates from runs i and j , fth being the threshold energy , and i(f(s ) <
g ) 
 being the indicator of the function f being less 
 than scalar g at point s. this indicates 
 greater stability and more physically meaningful  or at least 
 more easily interpreted  exploration .
the nonreapeatable free 
 energy basin predictions seen with unrestricted wtmetad do not correspond 
 to states that appear biologically relevant or well - populated at equilibrium 
 but rather appear to correspond to metastable traps found only because 
 of overzealous nonequilibrium kicks .
the mbmetad data proved more 
 useful for refining the choice of cvs and designing new simulations , 
 but because refining cvs is a system specific process rather than 
 a matter of mbmetad methodology , describing that process will be left 
 for future publication focused on the structure of g - actin . for our 
 purposes here , it is sufficient to demonstrate that the method enhances 
 sampling out of well - understood basins repeatably and without applying 
 risky nonequilibrium driving forces in less well understood regions 
 of cv space .
our final example concerns the clc - ec1 cl / h exchanger , a paradigmatic 
 protein 
 of the ubiquitous clc family . in this example
we will focus on investigating one 
 particular part of the molecular mechanism of chloride transfer : motion 
 of chloride anions in the channel when a key central acidic residue , 
 e148 , is protonated .
the previous examples showed that mbmetad is 
 an accurate and effective self - limiting mechanism when the cvs are 
 appropriate in a metastable basin .
however , it is often the case that 
 cvs will not be ideal even in the domain of interest , and with this 
 final example we intend to show that mbmetad continues to have both 
 desirable self - limiting behavior and dynamics with desirable hysteresis 
 in that case .
therefore , in this case our cvs are simply the displacements 
 along the membrane normal direction , z , from each 
 of two chlorides , one closer to the cytosolic side of the membrane , zcyt , and the other closer to the extracellular 
 side , zext , to the center of mass of a 
 central group of protein atoms that is chosen only as a stable reference 
 for the protein frame of reference and has no other intended physical 
 significance .
this leaves motion of e148 , water fluctuations , helix 
 motion , and other protein environment fluctuations as potential hsvs .
we aim here to compare the function of mbmetad with unrestricted 
 wtmetad in estimating pmfs of simultaneous translation of these chlorides 
 throughout the protein , with the mechanism presumed to involve five 
 key chloride states identified in previous structural studies .
first , when z is greater than 1 nm , the chloride is said to be extracellular . 

second , when z is near 0.6 nm , the chloride is expected 
 to be in its external binding site .
third , when z is approximately 0.5 nm , it is expected to be in its central 
 binding site . finally , when z is less than 1 , 
 this coordinate does not distinguish well between whether the chloride 
 is in its internal binding site or unbound on the cytosolic side of 
 the membrane .
the clc model system was set up as a dimer of 
 the wild - type clc - ec1 
 structure ( pdb : 1ots ) embedded into a lipid bilayer ( 163 pope ) and solvated with 11,000 
 tip3p waters in a 92  92  
 79  box under periodic boundary conditions .
protein 
 and lipids were modeled with the charmm27 force field ; long - ranged electrostatic interactions were treated with the particle 
 mesh ewald ( pme ) method , and the cutoff 
 distance for the short - ranged interactions from both the lennard - jones 
 and the real - space coulomb interaction was set to be 12  .
the 
 classical md simulation was performed using the gromacs md package patched with a customized version of plumed2 and a 2 fs time step .
the initial configuration 
 of the system was taken from our previous work , where the system was first equilibrated for 7 ns in the 
 npt ensemble with a temperature of 310 k and a pressure of 1 atm and 
 then further equilibrated for 10 ns in the nvt ensemble with a temperature 
 of 300 k. residues e113 in monomers a and b and d417 in monomer a 
 were protonated as suggested by previous calculations , and e148 in monomer a was also protonated to 
 allow chloride translations from the external to the central site . 

metadynamics was applied using the cvs described above for the 
 chloride ions in monomer a only , as previous experiments showed that each monomer carries out cl / h exchange independently .
monomer b was 
 still included in the system and simulated as normal , not subject 
 to any direct bias .
hills of initial height 0.2092 kj / mol and widths 
 0.035 nm in both dimensions were added each 1 ps . for wtmetad runs , 
 the well - tempered bias factor was 13 , whereas for transition - tempered 
 mbmetad runs the transition - tempered bias factor was 5 using transition 
 wells ( 0.45 , 1.17 ) and ( 1.75 , 0.12 ) .
mbmetad domains were targeted 
 to the transition barrier energy between those points plus 12 kj / mol 
 and were updated every 100 ps if the exterior bias level had changed 
 by 1 kj / mol . to restrict sampling to cases where the chlorides were 
 in the protein in the appropriate locations in the protein interior , 
 the coordinates of the extracellular and cytosolic - side chlorides 
 were both restrained .
the cytosolic - side chloride was restrained to 
 the region 3.24 < x <  0.93 ; 
 2.58 < y < 0.38 ; 0.09 < z < 2.00 , while the extracellular - side chloride was restrained 
 to the region 2.90 < x < 0.70 ; 2.94 
 < y < 0.21 ; 1.49 < z
the restraints were implemented as half - harmonic walls 
 with force constants of 200 kj/(mol  ) .
bias - based 
 free energy estimates from unrestricted wtmetad ( left 
 column ) and transition - tempered mbmetad ( right column ) over three 
 runs ( rows ) , each run until features stabilized qualitatively or appeared 
 to be irreparably incorrect . clockwise from top left ,
these pmf estimates 
 correspond to 250 , 300 , 300 , 300 , 350 , and 400 ns of simulation .
none 
 appear fully converged , yet the latter are more physically plausible 
 and the former show distinct signs of dynamics gone irrevocably astray 
 in the top and bottom runs ( see text for descriptions of the atomistic 
 details behind these unphysical features ) .
the free energy estimates given by three simulations using 
 each 
 method after various times , shown in figure 8 , do not correspond to converged potentials 
 of mean force .
they reflect the essential randomness of attempting 
 enhanced sampling while neglecting key hsvs and show that mbmetad 
 prevents the nonequilibrium driving from spoiling the dynamics as 
 it does so frequently in the unrestricted simulations . in the unrestricted 
 case ,
the metastable band structure seen in the mbmetad outside of 
 the starting basin is overshadowed by other features in two of the 
 three wtmetad replicates .
the top left energy profile shows 
 an unexpected and implausible 
 prediction of a sharp well where both chlorides are outside of the 
 protein that is much deeper than any of the expected binding sites , 
 while the bottom left shows a sharp line near zcyt = 1 that is not seen in any other simulations and also 
 does not correspond to any plausible physics : furthermore , the estimate 
 beyond this line appears to simply be a blob without the expected 
 basin structure for zext . in the simulation 
 corresponding to the top unrestricted wtmetad figure
, the trajectory 
 shows that the cystosolic chloride becomes stuck inside the protein 
 between helices that are normally firmly bound together . in the simulation 
 corresponding to the bottom unrestricted wtmetad figure ,
one of the 
 main helices unfolds and extends into the solvent : the chloride associates 
 with the helix , and then the nonequilibrium bias ratchets the helix 
 apart by pushing on the chloride , unfolding it turn by turn , before 
 the chloride dissociates and goes into the solvent , with effectively 
 permanent damage done . judging by the noticeably different energy 
 scale , it is evident 
 that mbmetad has expected self - limiting behavior even in this highly 
 challenging and poorly tuned test case , and judging by the fact that 
 the mbmetad runs do not exhibit the odd behavior seen in the others , 
 it also appears that self - limiting behavior improves robustness to 
 cv quality as predicted .
however , it is also clear that the mbmetad 
 runs in the right column of figure 8 do not fully agree among one another , showing that 
 it does not cover up core deficiencies in the cvs that can be seen 
 by comparing multiple independent replicates .
the disagreement 
 occurs at the level of which regions of interest 
 are selected on - the - fly , as well .
however , in the bottom right panel , the final region of interest 
 does not include the upper right quadrant of the allowed cv space . 
 at the time of termination ,
the simulation was sampling the upper 
 right region but not adding bias . given time ,
the histogram correction 
 would ensure that this region would properly register as part of the 
 cv space to bias .
still , as seen in the alanine dipeptide example , 
 it is better to use good cvs and a robust region definition than to 
 rely on long - time convergence of the domains .
this reflects the inefficiency 
 of any cv - based enhanced sampling in the presence of hsvs .
the bias - based pmf estimate of the primary trajectory at 200 ns ( left ) 
 is compared with estimates at 300 ns from direct continuation of that 
 trajectory ( middle ) and continuation from a new random equilibrated 
 initial conformation ( right ) .
the prominent differences between the 
 two show that sampling remains highly autocorrelated on a 100 ns time 
 scale even when biased , a hysteresis that properly indicates the presence 
 of hidden barriers .
next , mimicking the usual 
 process of checking whether one run might 
 have luckily converged to the correct result despite the gauntlet 
 of possible issues that evidently affected the other runs , we initialized 
 another mbmetad replicate of length 100 ns beginning with the bias 
 and stored histogram of the top run just before 200 ns to check if 
 the pmf estimate of this replicate would agree with the original copy 
 of the top run after 300 ns . as figure 9 makes evident , that is not the case .
the original 
 run and its branched replicate differ ; the former samples configurations 
 in which both chlorides are far from the center of the protein , whereas 
 the other samples configurations in which the cytosolic - side chloride 
 resides in the central pocket .
the new replicate primarily explores 
 the initial basin , whereas the original primarily explored the other 
 basins .
motion between the basins is not fully facilitated by the 
 bias because key hsvs are involved in the true barriers between them , 
 and the bias reports that faithfully .
mbmetad correctly reports the 
 presence of these hsvs even as it prevents the less physically interesting 
 and implausible dynamics seen in the unrestricted wtmetad simulations . 

the 
 new method mbmetad fundamentally eliminates one of the primary 
 limitations of metadynamics , that adding energy indefinitely causes 
 undesired and irreversible change in many sensitive dynamical systems , 
 by providing an effective and convenient self - limiting mechanism that 
 causes it to fill up to a flexibly defined free energy level and no 
 farther  requiring no a priori estimate of 
 barrier heights but rather estimating them on - the - fly as necessary . 

this allows for more focused study that should be especially practical 
 when using many cvs at once and incidentally demonstrates a solution 
 to the problem of boundary artifacts in metadynamics on domains of 
 any shape in any number of dimensions .
most importantly , it means 
 that using mbmetad makes designing good cvs for metadynamics simpler . 

if energy is no longer added in regions of cv space that are irrelevant 
 to transition mechanisms , the cvs no longer must be carefully tuned 
 there .
furthermore , having a controlled energy level makes setting 
 the 
 hill height and the tempering rates for wtmetad and ttmetad simpler . 

current practice requires users to choose their desired energy level 
 by balancing the level of tempering with the initial rate of biasing ; but the tempering must also be matched to only partially known cv 
 hysteresis time scales , and the biasing rate must also be matched 
 to the only partially known rate of dissipation in the system .
thus , 
 the parameters are overdetermined in a manner depending on unknowns 
 with opaque trade - offs ; using mbmetad to disentangle the choice of 
 final energy level from the choice of tempering rate and hill height 
 can make the choices more straightforward .
additional technical complexity 
 of the method for the implementer is therefore counterbalanced by 
 decreased intuitive complexity for actual application . by restricting 
 our changes to the hill function alone
, we ensure 
 compatibility with many other advances in metadynamics methodology 
 that can not be discussed in depth here .
for instance , the new procedure 
 is compatible in principle with untempered , well - tempered , and transition - tempered metadynamics , multiple walkers and bias exchange , driven 
 metadynamics , experimentally directed 
 ( ensemble - biased ) metadynamics , concurrent metadynamics , multiple time steps , and any choice of cv .
it is not compatible with adaptive gaussian 
 hills but is compatible with field - coordinate 
 metadynamics .
as noted earlier , it is 
 implemented in a public fork of the plumed2 package and is available with a short guide upon request .
moreover , there is little reason to think that this strategy applies 
 only to metadynamics .
all other adaptive enhanced sampling methods 
 such as the adaptive biasing force and orthogonal 
 space random walk approaches can similarly 
 suffer hysteresis effects related to adding reckless driving forces 
 when their parameters are set too aggressively .
it may be that adding 
 region - limiting mechanisms like the one presented here may allow for 
 more of those aggressive parameter choices to be used safely , with 
 substantial potential efficiency gains . in each case , one could use 
 the bias and an auxiliary histogram together to determine regions 
 of interest on the fly based on the simple definitions we propose 
 here , such as  everything less than 2 kt above the barrier 
 between these two states. finally , we find the physical 
 underpinnings of mbmetad unexpectedly 
 natural , and we hope they will inspire further thought in the field . 

the new hills can be understood as approximate green s functions 
 of diffusion inside the domain of interest plus an approximate quasistationary 
 distribution in the exterior of the domain .
mathematical work on the 
 foundations of accelerated dynamics is finding 
 deep theoretical power in using quasistationary distributions in lifting 
 the green s function of a coarse - grained markov process onto 
 a fine - grained configuration space .
therefore , the mbmetad hills appear 
 to emerge as approximate mixed - resolution green s functions , 
 connecting the worlds of markov state modeling , accelerated dynamics , 
 multiscale modeling , and cv - based adaptive enhanced sampling in a 
 surprising way .
 metadynamics 
 is an enhanced sampling method designed to flatten 
 free energy surfaces uniformly . 
 however , the highest - energy regions 
 are often irrelevant to study and dangerous to explore because systems 
 often change irreversibly in unforeseen ways in response to driving 
 forces in these regions , spoiling the sampling . 
 introducing an on - the - fly 
 domain restriction allows metadynamics to flatten only up to a specified 
 energy level and no further , improving efficiency and safety while 
 decreasing the pressure on practitioners to design collective variables 
 that are robust to otherwise irrelevant high energy driving . 
 this 
 paper describes a new method that achieves this using sequential on - the - fly 
 estimation of energy wells and redefinition of the metadynamics hill 
 shape , termed metabasin metadynamics . 
 the energy level may be defined a priori or relative to unknown barrier energies estimated 
 on - the - fly . altering 
 only the hill ensures that the method is compatible 
 with many other advances in metadynamics methodology . 
 the hill shape 
 has a natural interpretation in terms of multiscale dynamics , and 
 the computational overhead in simulation is minimal when studying 
 systems of any reasonable size , for instance proteins or other macromolecules . 
 
 three example applications show that the formula is accurate and robust 
 to complex dynamics , making metadynamics significantly more forgiving 
 with respect to cv quality and thus more feasible to apply to the 
 most challenging biomolecular systems . 
Introduction
Methods
Results
and Discussion
Conclusions