nucleoside triphosphates 
 ( ntp ) contain a nucleoside bound to three 
 phosphate groups .
the hydrolysis of ntps such as atp or gtp by proteins 
 is one of the most important reactions in biology .
it provides energy 
 to biochemical functions that are of great importance to life and 
 plays essential signaling functions .
for example , adenosine triphosphate 
 ( atp ) , one of the most important ntps , is often considered as a  molecular 
 unit of currency  in intracellular energy transfer and plays essential roles in biological processes 
 such as synthesis of proteins and membranes , movement of the cell , 
 cellular division , etc .
another important ntp molecule , 
 guanosine triphosphate ( gtp ) , is vital to cell signaling , cell motility , 
 protein synthesis and translocation , and more . due to the importance 
 of ntp hydrolysis
, great experimental effort has been devoted to study 
 the hydrolysis mechanism for many years .
these efforts include but 
 are not limited to the structural determination of nucleotide - bound 
 proteins through x - ray crystallography and cryo - electron microscopy , as well as the identification of essential residues 
 through mutagenesis and biochemical assays . despite the success of such experiments , the specific molecular 
 mechanism(s ) of ntp hydrolysis often remains ambiguous because experiments 
 can not directly observe chemical reactions in proteins . as a result , 
 theoretical studies , conducted through modeling and simulation , are 
 highly desirable to reveal the hydrolysis mechanism at an atomistic 
 level .
even though molecular dynamics ( md ) simulations have 
 benefited 
 enormously from the rise of computational capabilities over the last 
 few decades , the accurate computational 
 study of ntp hydrolysis remains very challenging for primarily two 
 reasons .
first , the nature of the hydrolysis reaction requires the 
 chemical bonding to be modeled in a reactive manner ( i.e. , bonds will 
 be broken and new ones will be formed ) .
quantum mechanics / molecular 
 mechanics ( qm / mm ) simulation is one possible solution to this problem , 
 in which the ntp molecules and nearby protein region undergoing and 
 directly influencing the hydrolysis reaction are modeled with electronic 
 structure theory , and the other parts of the system are modeled with 
 classical molecular mechanics . however , the phosphate groups , catalytic 
 metals , and their nearby amino acids and water molecules are all relevant 
 to the hydrolysis reaction , usually making the reactive qm region 
 larger than 150 atoms . simulating
a system of this size with electronic 
 structure theory is extremely computationally expensive , even at the 
 level of density functional theory ( dft ) . to overcome this challenge ,
warshel and co - workers have proposed the empirical valence 
 bond ( evb ) approach , which uses a diabatic representation 
 to describe the reactant , intermediate , and product states of the 
 reaction and empirical interaction , as one way of studying reactions 
 in complex systems .
the multiscale 
 reactive molecular dynamics ( ms - rmd ) approach , where the reactive 
 potential energy surface is defined by an  on the fly  
 dynamically evolving linear combination of multiple diabatic potentials , 
 and the terms of the interaction matrix defining them are variationally 
 determined ( at least in part ) from qm / mm or ab initio md forces , could in principle be applied to ntp reactions , as has 
 been done to simulate proton transport in proteins ( see , e.g. , ref ( 17 ) for a recent example ) . 

the ms - rmd reactive potential energy surface enables the formation 
 and cleavage of molecular bonds in a deterministic classical md simulation . 
 as another alternative , the self - consistent charge density functional 
 tight binding ( scc - dftb ) might be utilized . 

this approach is a semiempirical qm algorithm that is derived from 
 a taylor expansion of the dft kohn 
sham energy for charge density 
 fluctuations in which the hamiltonian matrix elements are evaluated 
 using a minimal basis set of pseudoatomic orbitals within a two - center 
 approximation .
another possible approach 
 is to utilize force matching along a reaction pathway determined from 
 less expensive semiempirical calculations such as scc - dftb and to 
 then correct them to higher level qm calculations .
in general , all of the aforementioned methods are implemented 
 in order to greatly increase the amount of statistical sampling in 
 the simulation by making evaluations of the forces in the qm region 
 of the hydrolysis reaction faster , ideally ( but not commonly ) with 
 a systematic and small trade - off in accuracy .
the second challenge 
 is that the free energy barrier of ntp hydrolysis 
 is usually on the order of tens of kcal / mol , which makes conventional 
 md and qm / mm simulations of no practical use in calculating a potential 
 of mean force ( pmf ) because the transition state region will barely 
 be sampled during a feasible simulation time . in a recent qm / mm study 
 of atp hydrolysis in f1-atpase ,
shigehiko 
 et al . optimized the geometry of the hydrolysis compounds along the 
 reaction path and constructed an approximate free energy surface with 
 vibrational frequencies from qm calculations .
this approach only provides 
 free energy profiles for the stationary points ( potential energy minima 
 and transition states ) instead of the entire reaction path , and , more 
 importantly , the lack of accounting for configurational entropy from 
 nonlinear fluctuations may cause significant errors in the free energy 
 surface for enzymatic reactions . therefore , 
 enhanced sampling methods have been employed in the studies of ntp 
 hydrolysis to address this high energy barrier issue .
used umbrella sampling and 
 the weighted histogram analysis in their scc - dftb / mm 
 simulations to encourage the attack of lytic water to the gamma phosphate 
 and the proton transfers between the lytic water and amino acid .
more 
 recently , metadynamics ( metad ) has been adopted in ntp hydrolysis 
 simulations .
metad accelerates the simulation by discouraging 
 it from revisiting configurations that have already been sampled through 
 constructing a history - dependent bias potential v(s ) to help the system to overcome large free energy 
 barriers .
an appealing feature of metad is its ability to explore 
 reactions without an a priori guess of the detailed 
 reaction path , in contrast to umbrella 
 sampling where one has to prepare the geometry for each window according 
 to a hypothetical reaction path . in spite of the past application 
 of metad in assisting the pmf 
 calculation of ntp hydrolysis ,
the practical and robust convergence 
 of an acceptably accurate and therefore computationally expensive 
 qm / mm simulation with metad remains extremely difficult .
other than 
 the sampling limitations of the qm / mm simulation discussed earlier , 
 the complexity of ntp hydrolysis may demand more collective variables 
 ( cvs ) than those solely associated breaking the bond between the beta 
 phosphate and the gamma oxygen , e.g. , the proton transfer between 
 the lytic water and/or the amino acid through a water mediated  proton 
 relay  .
additionally , different 
 sets of cvs are sometimes required to properly sample the forward 
 and backward reactions . as a result , 

the quality of the convergence in these early qm / mm metad simulations 
 could not be verified as rigorously as the theory suggests , e.g. , 
 by checking for the diffusive motion of the cvs and for agreement 
 between multiple replicas .
it is noteworthy that these early 
 qm / mm metad simulations have 
 either utilized the original nontempered metad , where the same gaussian is added to the hamiltonian along 
 the simulation , or the well - tempered metad ( wtmetad ) , where the gaussian is scaled with respect to the current 
 local bias potential . in nontempered metad
, the ever - growing history 
 dependent bias potential does not converge to the pmf modulo a constant 
 but rather fluctuates around the desired final state .
a possible mitigation of this behavior is to take the time 
 average of the bias potential after observing the diffusive motion 
 of the system in the cv space .
however , as pointed out earlier , the 
 identification of diffusive motion can be subjective and itself is 
 very difficult to achieve due to the limited sampling of ntp hydrolysis 
 in expensive qm / mm simulation , thus making an extensive enough metad 
 simulation to check for truly diffusive behavior is presently nearly 
 impossible .
wtmetad overcomes the convergence issue but presents an obvious trade - off between the exploration 
 of the cv space in the early stages of the simulation and the asymptotic 
 convergence rate .
for example , wtmetad that decreases 
 the gaussian height slowly builds up a fast - growing but noisy bias 
 potential , hindering the convergence .
in contrast , a fast decrease 
 of the gaussian heights generates a smooth bias potential but prevents 
 it from escaping a minimum within a feasible qm / mm simulation time , 
 especially when the barrier in the cv space is very large , which is 
 usually the case in ntp ( e.g. , atp ) hydrolysis .
as demonstrated 
 in our recent papers , the newly developed transition - tempered 
 metad ( ttmetad ) has been shown to converge 
 rapidly and asymptotically 
 without sacrificing the exploration of the cv space in the early stages 
 of simulations , in contrast to current other convergent metad methods . 

this aspect of ttmetad obviously makes it very appealing to study 
 ntp hydrolysis in computationally expensive qm / mm simulations , as 
 the sampling of the latter will always be limited due to the high 
 cost of the qm force evaluations . additionally , as pointed out by 
 mcgrath et al .
, as well as in our own 
 experience , in a longer metad simulation the bias potential may introduce 
 irreversible changes in the system that dramatically delays the convergence 
 or pushes the system into unstable regions and may even crash the 
 qm / mm simulation .
therefore , it is also of interest to employ a self - limiting 
 bias potential to avoid these issues and to properly focus the qm / mm 
 free energy sample , by using a method such as our recently developed 
 metabasin metad ( mbmetad ) . the 
 purpose of this paper is to demonstrate , through the specific 
 example of atp hydrolysis in actin , that a combination of ttmetad 
 and mbmetad with qm / mm md simulation can provide a systematic and 
 rigorous estimate of the free energy pathway of ntp hydrolysis in 
 proteins .
we study here the atp hydrolysis in monomeric and filamentous 
 actin and compare the features of the two reactions , building on our 
 earlier work with this more powerful 
 and convergent combinations of simulation methods .
atp hydrolysis 
 in actin plays an important role in the actin - based cytoskeleton .
for example , atp - bound monomeric globular actin ( g - actin ) is added 
 to the ( + ) end of an actin filament , undergoing conformational change 
 to a more flattened state in filamentous actin ( f - actin ) , which then 
 catalyzes atp hydrolysis , then releases the hydrolyzed phosphate group , 
 and finally dissociates at the (  ) end , now as an adp - bound 
 g - actin .
the high selectivity for hydrolyzing atp in f - actin , as opposed 
 to in g - actin , is important for modulating the filament s physical 
 properties and the binding affinity of various actin - binding proteins . 

the most recent experimental studies have reported 
 the atp hydrolysis rate to be 0.3  0.1 s and 7  10 s for the 
 reaction in f - actin and g - actin , respectively .
assuming the same arrhenius 
 prefactor in both forms of actin , this dramatic speedup suggests a 
 6.6  0.7 kcal / mol hydrolysis free energy barrier difference 
 between f- and g - actin .
mccullagh et al . carried out a qm / mm md computational study to reveal the origins 
 of the different rates of atp hydrolysis in the different actin states . 

nontempered metad was employed in that work , and the simulations showed 
 an 8 kcal / mol hydrolysis barrier difference between f- and g - actin , 
 in good agreement with the experimental estimate .
the water structure 
 of the  proton relay  between the lytic water in the 
 hydrolysis site and asp154 was found to be to one difference between 
 f - actin and g - actin . in this paper
, we revisit these same actin 
 systems and compare 
 the performance of ttmetad simulation against wtmetad and nontempered 
 metad simulation .
our simulation follows the setup of mccullagh et 
 al . , except for using different metad 
 approaches .
in particular , we show that reaching convergence requires 
 a comprehensive representation of the overall hydrolysis process by 
 also adding in an additional cv that biases proton transfer from the 
 nucleophilic water molecule .
this additional concurrent metadynamics 
 benefits from the self - limiting feature of mbmetad and dramatically 
 accelerates the convergence of the free energy surface .
most importantly , 
 our simulations provide a rigorous analysis of the convergence of 
 the bias potential , which has so far largely been overlooked in qm / mm 
 simulations of ntp hydrolysis .
the focus and conclusions of this paper 
 reveal the computational advantages of utilizing ttmetad and mbmetad 
 with qm / mm
future work 
 will focus on utilizing this approach to then study the hydrolysis 
 mechanism itself in key proteins .
metad utilizes an 
 external history dependent bias potential , which has the form of a 
 summation of gaussian functions centered at previously visited points 
 in a preselected cv space .
the bias potential helps the chemical system 
 to overcome large free energy barriers by discouraging it from revisiting 
 configurations that have already been sampled . in the first implementation 
 of metad , the bias potential ( vg ) evolved 
 according to the equation1with w called the bias energy 
 addition rate ,
usually expressed as the height of individual gaussian 
 functions added , w0 , and divided by a 
 deposition stride . the function si(r ) represents the i cv for metad simulations , i.e. , distance , angle , dihedral , 
 atomic coordinates , radius of gyration , local density , etc . the term f(si(r ) ) is defined as2 in which i is the width of the gaussian for the i cv .
this approach is usually referred to 
 as nontempered metad because 
 the bias energy addition rate w does not depend on 
 the progress of the simulation .
the bias potential from nontempered 
 metad does not converge asymptotically , and a practical solution is to take the time average of the bias 
 potential , after the time ( t0 ) when the 
 motion of the system on the cv space becomes diffusive .
the pmf estimated from nontempered metad takes 
 the form3where vg(t ) is the instantaneous 
 bias potential . to address 
 the convergence issue
, wtmetad tempers the gaussian heights exponentially 
 with respect to the local bias potential .
the height of the gaussian 
 function becomes time - dependent4 in which t is a parameter 
 ( with the unit of temperature ) chosen prior to the simulation .
dama 
 et al . have provided a rigorous mathematical proof that wtmetad converges 
 asymptotically to a linearly scaled inverse of the underlying pmf . on the other hand , ttmetad smoothly converges 
 the bias potential 
 asymptotically like wtmetad , but it keeps track of the region that 
 the bias potential has explored and delays the tempering of the gaussian 
 height until all the basin points are included in this region .
because 
 the basins are relatively filled prior to the tempering , ttmetad can 
 be used with a quickly converging tempering ( small t ) without a concern for being trapped in the first nearby 
 free energy minimum .
the height of the gaussian in ttmetad is tempered 
 as5where v*(s , t ) is 
 the minimal bias on the maximally biased path among all the 
 continuous paths s( ) that connects all of 
 a set of preselected basin points in the cv space .
instead , one needs to estimate the position of the basin 
 points in the cv space .
this estimation is straightforward when it 
 is defined by the problem of interest , i.e. , in this case the atp 
 state and adp + pi state for the atp hydrolysis reaction .
generally , 
 ttmetad has been proven to be robust to the accuracy of the basin 
 points as long as they are separated by a major barrier . in a conventional 
 metad simulation
, 
 the bias energy is constructed to flatten the underlying free energy 
 surface globally on the cv space . however , in many - dimensional systems , 
 the system may be pushed out from the low energy regions that are 
 of interest to some of the higher energy regions that are not physical 
 or as relevant .
worse , entering these high - energy regions can couple 
 with large , effectively irreversible conformation changes , introducing 
 large hysteresis and limiting reproducibility in the estimates of 
 the pmf . therefore , it is highly desirable to use a self - limiting 
 biasing method that only fills the underlying free energy surface 
 up to a specified level and then converges that bias potential to 
 estimate the pmf .
this forces a more physically relevant exploration 
 and greater reproducibility of the pmf estimates .
mbmetad achieves 
 this self - limiting behavior for the bias potential by only flattening 
 a finite domain ( d ) of the whole cv space .
mbmetad 
 therefore utilizes a gaussian with two parts6 in which fint(si(r ) ) is the 
 part to flatten the interior of the domain , and fext(si(r ) ) raises the bias level of the exterior to match the bias 
 level of the domain boundary to counteract the force from fint(si(r ) ) that slopes into the exterior of the 
 domain . the term fint(si(r ) ) is rescaled from 
 the regular gaussian ( named f0(si(r ) ) ) following 
 a relative to the multiplicative mcgovern - de pablo rule7with i being the boundary - normalized 
 integral of f0(si(r ) ) over d , and f is a scaling function .
the function f(x ) evaluates as x for x  1 and x + 0.5*(1x ) for x
< 1 , so that these 
 hills are identical to the mcgovern - de pablo form inside the domain 
 but taper more smoothly to zero outside of the domain . at the same 
 time
, mbmetad constructs the exterior bias update , fext(si(r ) ) , as the function that would force the sampling into 
 the domain near the boundary exactly enough to cancel out how the 
 interior bias , fint(si(r ) ) , would push the 
 sampling out near the domain boundary if sampling inside the domain 
 were perfectly flat ( or were the desired tempered distribution in 
 the case of wtmetad ) .
a detailed description of mbmetad and additional 
 applications can be found in ref ( 37 ) .
the structures 
 of both g- and f - actin 
 were generated from equilibrated classical molecular dynamics simulations .
the crystal structure ( pdb code 1nwk ) was used for the 
 g - actin simulation and the oda model ( pdb code 2zwh ) after equilibration 
 of a 13-mer periodic filament was used for f - actin simulation .
the 
 main change during the actin polymerization is a flattening of the 
 subdomain 2134 dihedral , which is 0.55 
 degree in f - actin and 27.09 degree in g - actin . to prevent 
 f - actin from relaxing to g - actin during the metad simulation , a harmonic 
 bias force with a spring constant of 228.8
kcal / mol was applied to 
 restrain this dihedral based on a coarse - grained mapping of f - actin . 

further details of the structures in this simulation are provided 
 in mccullagh et al . and its supporting 
 information
, including methyl triphosphate , 
 bound magnesium , surrounding waters , and ten amino acids ( asp11 , gly13 , 
 ser14 , gly15 , lys18 , gln137 , asp154 , gly156 , asp157 , and his161 ) . 
 the atoms involved in coordination with the catalytic magnesium , i.e. , 
 one beta oxygen , one gamma oxygen , and four water molecules , are included 
 in this qm region .
sequential amino acids backbone atoms were included 
 and hydrogen capped at the terminal nitrogen and carbon .
nonsequential 
 amino acids were truncated at the c or c carbon and capped with hydrogens .
the 
 pbe / tzv2p level of density functional theory was employed for the qm region simulations , and the magnesium ion 
 was treated with gth pseudopotentials .
the charmm27 
 force field was used to model the mm region , 
 and periodic boundary conditions were applied in the nvt ensemble 
 with a smoothed particle mesh ewald treatment of the long - range electrostatics .
the temperature was maintained at 310 k with 
 a nos  hoover thermostat , and the integration time step was set to 0.5 fs , totaling over 1.5 
 ns of simulation time .
the qm / mm metad simulations were performed 
 with the cp2k software coupled to the 
 enhanced sampling package plumed2 .
a set of appropriate cvs ( denoted 
 as  s  ) representing the chemical reaction 
 of interest must be chosen prior to metad simulation .
as described 
 earlier , metad sequentially builds up a history dependent bias potential 
 of these variables to accelerate their dynamics . as suggested by mccullagh 
 et al .
, our ttmetad and wtmetad simulations 
 utilize two cvs in the atp hydrolysis reaction ( cf .
one cv is the coordination 
 number between the terminal ( gamma ) phosphorus ( p ) and bridging ( beta ) oxygen ( o ) , shown in figure 1 .
this coordination 
 number describes the forming and cleavage of the bond between p and o and allows the recombination 
 of p with any of the o atoms . 

the other cv , also seen from figure 1 , is the coordination number between p and both gamma oxygen ( o ) and qm water oxygens 
 ( ow ) .
therefore , all the qm water molecules ( about 30 ) 
 can possibly be the lytic water that attacks p and 
 triggers the hydrolysis reaction , in contrast to only biasing a particular 
 water molecule in the qm region to be the lytic water .
we anticipate this species - agnostic approach 
 to be important to prevent unphysical artifacts from breaking chemical 
 symmetries while modeling this reaction .
an illustration of the 
 atoms involved in the two collective variables 
 employed in the ttmetad simulations of atp hydrolysis in actin .
it should be noted that the cvs 
 in our simulation prefer neither 
 an  associative  mechanism nor a  dissociative  mechanism for the atp hydrolysis .
the former suggests the lytic water molecule 
 forms a bond with p before the bond between p and o breaks .
the latter proposes 
 the p and o bond breaks before 
 the nucleophilic attack of lytic water .
importantly , an additional 
 cv that biases the hydrogen transfer from the lytic water was also 
 used in our qm / mm ttmetad simulation with a self - limiting metad ( mbmetad ) 
 to enhance convergence of the pmf .
the parameter t ( in eq 4 ) 
 that tunes the speed of decreasing gaussian 
 height significantly affects wtmetad convergence rate .
it is suggested in the original wtmetad paper that t should be chosen according to the free energy barrier , 
 as t + t should be of the 
 order of magnitude of the barrier height .
therefore , t was chosen to be 12090 k ( corresponding 
 to 24.0 kcal / mol ) to maximize the efficiency of wtmetad , based on 
 the previously determined hydrolysis barrier height ( 30 and 22 kcal / mol 
 for atp hydrolysis in g- and f - actin , respectively ) .
however , it is important to note that this level of accuracy 
 of the estimates of the energy barrier is almost impossible for every 
 npt hydrolysis reaction , especially if one is interested in a reaction 
 that has not been studied before .
we note that ttmetad is more robust 
 to the choice of t as compared to wtmetad , and it is recommended to temper more aggressively ( with a smaller 
 t value ) , since the tempering is not applied 
 until the basins are already relatively filled .
we therefore set t to be 1240 k ( corresponding to 2.5 kcal / mol ) for the ttmetad 
 simulations .
additionally , one needs to provide the estimates of the 
 basin positions , which are easy to do for atp hydrolysis ( one basin 
 is the atp region and the other is the adp + pi region ) .
the same 
 heights of gaussian are used for both wtmetad and ttmetad simulations 
 as in the previous nontempered metad simulations .
we carried out three sets of qm / mm 
 simulations to compare the performance between wtmetad and ttmetad 
 in the early stages of the simulation : f - actin starting from atp , 
 f - actin starting from adp + pi , and g - actin starting from atp . for 
 hydrolysis in f - actin , wtmetad and ttmetad simulations were initiated 
 separately from atp and adp basins . to illustrate
the data , the bias 
 potentials from the atp and adp replicas were summed and plotted together 
 into one plot ( top row in figure 2 ) , following the practice of multiple walker metad .
this figure demonstrates the efficiency of ttmetad 
 in exploring the cv space at the early stages of the simulation .
it 
 is seen that after 50 ps , both atp and adp + pi basins of the underlying 
 free energy surface are much better filled in ttmetad simulations 
 than in wtmetad simulations . in the ttmetad simulations , both atp 
 and adp replicas have escaped from their origin basin and started 
 to explore the transition state region , in contrast to them being 
 trapped in its origin basin in wtmetad simulations .
the wtmetad contour 
 lines are much smoother than those in the ttmetad plot , reflecting 
 the fact that the wtmetad is using smaller gaussians before overcoming 
 the barrier .
the maximum value of the bias energy , which can be identified 
 from the background color in each figure , is about 25 kcal / mol in 
 ttmetad , significantly higher than the value of 18 kcal / mol in wtmetad . 
 according to the tempering rules of wtmetad ( eq 4 )
, a continued wtmetad simulation would use 
 as small as 47% of the original gaussian hill height , in contrast 
 to ttmetad that still uses a full gaussian since the predefined basins 
 have not been connected by the bias potential yet .
wtmetad will eventually 
 converge with longer simulation times ; however , due to tempering , 
 it will take even longer to add the same amount of bias energy to 
 the system than before , and this makes it relatively inefficient . 

the g - actin atp simulations show a very similar story to the f - actin 
 simulations ( bottom row in figure 2 ) .
inverse of the bias potential for atp hydrolysis in f - actin 
 ( top 
 row ) and g - actin ( bottom row ) . for f - actin ,
the bias potential is 
 the summation of atp and adp+pi replicas ( 50 ps each ) . for g - actin , 

bias potentials 
 from wtmetad are shown in the left column , and bias potentials from 
 ttmetad are shown in the right column .
the comparison 
 between ttmetad and wtmetad was clear , 
 so we halted the wtmetad simulations at 50 ps .
it is important to 
 emphasize that wtmetad was set up in a way to maximize its efficiency 
 by taking advantage of an accurate prior knowledge of the hydrolysis 
 barrier height , but ttmetad s 
 advantages seen early in the simulation were compelling .
indeed , the 
 early exploration efficiency of ttmetad appears to be of considerable 
 value to the computationally demanding qm / mm ntp hydrolysis simulations .
we 
 next continued all three sets of ttmetad simulations ( f - actin starting 
 from atp , f - actin starting from adp + pi , and g - actin starting from 
 atp ) beyond 50 ps to achieve better pmf convergence .
we also replicated 
 each set , so that the hydrolysis in g - actin and in f - actin was both 
 modeled by 4 replicas each .
after reading in the ttmetad bias potential 
 as shown in figure 2 , the replicated trajectories were restarted from the configuration 
 at 50 ps with a random velocity sampled from a boltzmann distribution 
 at 310 k. the reported pmf estimates are then averaged over those 
 replicas , because this ensemble - average approach has previously proven 
 to converge the metad simulations faster in both wall times and cpu 
 times .
as the ttmetad simulation 
 goes on , the trajectories show diffusive motion between atp and adp 
 + [ h2o  po3 ] .
however , 
 a complete hydrolysis reaction requires one of the protons from the 
 lytic water in [ h2o  po3 ] to transfer to o to form [ h2po4 ] ( see figure 3 ) ; but , without an additional bias , [ h2o  po3 ] ( seen in figure 3 ) remains a metastable state , and the transition between 
 [ h2o  po3 ] and [ h2po4 ] does not happen spontaneously .
the [ h2o  po3 ] and [ h2po4 ] correspond to a very similar value on the 2d cv space ( with only 
 slight differences due to surrounding water structure ) , and , as a 
 result , the 2-cv biasing force from the previously described ttmetad 
 simulations did not accelerate this proton transfer process .
the neglect 
 of an important slow degree of freedom ( i.e. , rare event ) in metad 
 simulations usually leads to significant hysteresis in the pmf estimates , and we observed this in our simulation . taking 
 the atp hydrolysis in g - actin as an illustrative example , the longer 
 simulation shows the trajectories spend a relatively longer time between 
 the [ h2o  po3 ] and 
 the transition state region , leading to an unphysical minimum near 
 that region in our pmf estimates ( right panel in figure 4 ) .
an illustration of the 
 collective variables employed in the mbmetad 
 simulations to accelerate the transition between [ h2o  po3 ] ( bound water proton coordination number 
 is 2 , shown as purple ) and [ h2po4 ] ( bound water proton coordination number is 1 , not shown ) .
bias potential from simulations of atp hydrolysis 
 in g - actin . left 
 panel : the inverse of the bias potential from averaging over 4 ttmetad 
 replicas with 2 cvs ( 180 ps for each replica ) .
right panel : the extended 
 ( 50 ps ) simulation from the left panel .
an unphysical hysteresis 
 builds up in the pmf ( 0.1 , 3.8 ) because of the neglect 
 of a third key slow degree of freedom ( proton transfer of the lytic 
 water ) .
the unit of energy is kcal / mol . in general , metad simulations with more than 2 instantaneous 
 cvs 

one solution to this higher - dimensional 
 problem is to simulate a series of metad simulations with different 
 concurrent cvs and let their bias energy and configuration exchange 
 between one another .
however , the application 
 of this approach to the atp hydrolysis system seems less than optimal , 
 as the number of replicas needed to build  the ladder of replicas  
 could be large and qm / mm simulation is costly
nevertheless , 
 we decided to explore the idea of concurrent metad and use it to explore just a third cv to accelerate 
 the hydrogen transfer of the lytic water : the coordination number 
 between lytic water oxygen ( ol ) and lytic water hydrogen 
 ( hl ) .
this cv was chosen to bias the transition between 
 [ h2o  po3 ] ( proton 
 coordination number of the bound water is 2 ) and [ h2po4 ] ( proton coordination number of the bound 
 water is 1 ) .
an illustration of this cv can be found in figure 3 . as shown in other research , 
 as well as in our own simulations ,
the lytic water , instead of an 
 already dissociated hydroxide , nucleophilically attacks 
 the backside of p and initiates the hydrolysis reactions . 

therefore , biasing the lytic water from the beginning of the simulation 
 has a risk of introducing artifacts into the simulation , i.e. , forcing 
 hydroxide to be the attacking group , if the metad drives autoionization 
 too aggressively
. additionally , the lytic water could be any water 
 molecule adjacent to the p , and biasing all of them 
 to the point of autoionization is clearly undesirable .
therefore , 
 we utilized mbmetad to introduce a concurrent and self - limiting 
 bias potential that features the following properties : 1 ) it is only 
 in play after the formation of the [ h2o  po3 ] intermediate ; 2 ) it applies only to the lytic 
 water that forms the [ h2o  po3 ] intermediate ; and 3 ) the maximum bias value applied to the cv is 
 less than 6 kcal / mol ( self - limiting ) .
the value of 6 kcal / mol was 
 chosen as an upper - limit so that the ol  hl bond would not break unless there is a good hl accepting 
 group close by , in contrast to conventional metad where the ol  hl bond would be driven to break regardless 
 of the local environment . as a result of this mbmetad bias ,
the weakened 
 ol  hl bond had a much shorter lifetime 
 than in the 2-cv ttmetad simulations , and the transfer between [ h2o  po3 ] and [ h2po4 ] became diffusive .
then , depending 
 on the chemical environment , the proton transfers from lytic water 
 oxygen to gamma oxygen followed different pathways .
ttmetad 
 combined with concurrent mbmetad is able to converge the 
 simulation well , and the final pmf estimates from averaging over four 
 replicas are depicted in figure 5 .
the minimum free energy path was then computed with 
 the string method and is plotted in figure 6 .
the difference 
 in activation free energy of atp hydrolysis between g - actin and f - actin 
 is 7.1 kcal / mol , aligning well with the experimentally suggested difference 
 of 6.6  0.7 kcal / mol .
this agreement is also closer than the 
 previous result of 8 kcal / mol obtained from nontempered metad simulations .
the overall magnitude of the hydrolysis barrier 
 is almost certainly dominated by electrostatic effects in the active 
 site as well as the energy required to break the chemical bonds plus 
 differences in the efficiency of the proton relay from the lytic water 
 in f - actin versus g - actin .
atp hydrolysis 
 in f- ( left ) and g - actin ( right ) pmf from ttmetad 
 + concurrent mbmetad simulations as described in the main text ( averaged 
 over 4 replicas ) .
minimum free energy path from the 2d free energy surfaces ( figure 5 ) for atp hydrolysis 
 in f- and g - actin .
the reaction coordinate represents the hydrolysis 
 reaction progress along the minimum free energy path , i.e.  0  
 being atp and  1  being adp + pi .
the asymptotic convergence 
 for tempered metadynamics ( wtmetad and ttmetad ) can not in principle 
 be reached with finite sampling .
therefore , 
 it is important to set up practical convergence criteria that can 
 be verified for expensive qm / mm simulations such as the present ones 
 of atp hydrolysis . as suggested in a number of papers , the diffusive motions of the cvs and the differences between the 
 replicas can be employed as convergence criteria in metad simulations . 

conversely , the evolution of the bias potential is not used because 
 of a concern for false convergence , where the gaussian height could 
 have been tempered so little that the bias potential barely changes 
 with longer simulation , but the calculation is in fact not converged . to verify the diffusive motion of the system in cv space , we restarted 
 trajectories from the atp configuration with velocities randomly sampled 
 from a boltzmann distribution ( at 310 k ) .
 
 pmf ( figure 5 ) was 
 read in , and the gaussians kept being added following the ttmetad 
 tempering rule .
the concurrent and self - limiting bias accelerating 
 the transfer between [ h2o  po3 ] and [ h2po4 ] was applied 
 as described before .
if the previous simulation had indeed reached 
 a reasonable degree of convergence , one would expect the system to 
 escape from the atp basin , cross the transition state region , form 
 adp + [ h2o  po3 ] , and 
 eventually reach adp + [ h2po4 ] relatively quickly , i.e. , the whole process should be diffusive 
 without being trapped in any region .
taking f - actin as an example , 
 we tested eight such trajectories , and all of them show satisfactory 
 diffusive motion ( i.e. , by finishing the hydrolysis process within 
 15 ps , which includes constructing the self - limiting bias for the 
 hydrogen transfer cv ) .
an example trajectory of atp hydrolysis in 
 f - actin that demonstrates the proton transfer is shown in figure 7 .
the evolution of 
 the cvs for another testing trajectory of hydrolysis in f - actin is 
 depicted in the supporting information ( figure 
 s1 ) .
we also compared the differences 
 between the bias potentials from 
 each replica because the robust convergence can not be confirmed until 
 they agree reasonably well with every other .
as mentioned earlier , 
 we simulated four replicated simulations for atp hydrolysis in both 
 f- and g - actin .
figure 8 depicts the minimum free energy paths ( the atp region is shifted 
 to zero ) from the four replicas of atp hydrolysis in g - actin .
all 
 four minimum free energy paths show good agreements in the hydrolysis 
 barrier .
however , differences in the adp + pi region indicate that 
 the convergence quality is not as ideal as in the atp transition state 
 region .
nevertheless , in actin the forward atp hydrolysis reaction 
 ( forming adp + pi ) is much more important than its synthesis ( reverse 
 reaction ) , which is a very rare event
. minimum free energy path 
 from 4 replicas of atp hydrolysis in g - actin . 

the reaction coordinate represents the hydrolysis reaction progress 
 along the minimum free energy path , i.e.  0  being atp 
 and  1  being adp + pi .
it is also of interest to compare the converged ttmetad and 
 concurrent 
 mbmetad simulation with the previous nontempered metad simulation 
 from mccullagh et al .
the pmf estimates 
 and the minimum free energy path from the nontempered simulation are 
 replotted in the supporting information ( figures s2 and s3 ) .
even though there are differences in the pmf 
 estimates between the two studies , the key feature , the barrier height 
 in the minimum free energy paths from both simulations , agrees well . 

as has been discussed in the introduction , 
 a qm / mm study of ntp hydrolysis is extremely expensive and complicated ; 
 therefore , mccullagh et al . had chosen 
 to employ nontempered metad to ensure the exploration of the cv space 
 rather than to aim for a more meaningful pmf convergence with a method 
 like wtmetad .
also , there was a lack of averaging of the nontempered 
 metad bias after observing the diffusive motion ; hence the convergence 
 of the nontempered metad pmf estimate was not ideal .
however , as shown 
 here with the recently developed ttmetad and concurrent mbmetad , qm / mm 
 simulations are now able to provide a rigorous and converged pmf estimate 
 for this system and presumably others similar to it .
protein mediated npt hydrolysis such as the actin atp hydrolysis 
 studied in this work remains one of the most important and challenging 
 biochemical reactions in computational biophysics .
the size and complexity 
 of the protein environment adjacent to the hydrolysis site ( water 
 structure , catalytic metal , and amino acid residues ) , as well as the 
 chemical complexity of the hydrolysis reaction itself , limits the 
 amount of sampling that one can achieve in a computationally demanding 
 qm / mm simulation .
therefore , an enhanced sampling method featuring 
 high - quality pmf estimates at the earliest possible stage of the simulation 
 is highly desirable . as shown in this work ,
ttmetad is well suited 
 to this task in that it explores the cv space with full gaussians 
 until the basins are roughly filled .
. 
 the ttmetad qm / mm simulations also exhibit diffusive motion between 
 atp and adp + [ h2opo3 ] .
we further introduced an additional bias potential 
 to accelerate 
 the transition between [ h2o  po3 ] and [ h2po4 ] . to perturb the hydrolysis dynamics to the least effective extent
it 
 is important to emphasize that our metad simulations were not carried 
 out with 3 simultaneous cvs , where one must compute the pmf with respect 
 to any two cvs at the fixed value of the third cv .
the third cv ( biasing 
 proton transfer from lytic water ) was only introduced to accelerate 
 one slow degree of freedom that prevented the ttmetad with other two 
 hydrolysis cvs from converging .
we therefore consider ttmetad + concurrent 
 mbmetad to be a much more efficient and physically motivated approach 
 than metad simulations with 3 simultaneous cvs .
the reasons are 1 ) 
 the larger number of simultaneous cvs that metad uses , the longer 
 the simulation takes to converge ; 2 ) biasing the lytic water molecules 
 at the outset of the simulation can lead to cleaving the o  h 
 bond before it nucleophilically attacks atp ; and 3 ) identifying in 
 advance a particular water molecule to be the lytic water from all 
 of the qm water molecules is nontrivial , so it is optimal to allow 
 one to initiate the reaction before adding the mbmetad bias to it .
we have also shown that running replicated tempered metad simulations 
 is advantageous , since not only do they provide confirmation of convergence , 
 but they can also reduce the simulation wall clock times , thus allowing 
 them to converge faster . furthermore , a single tempered metad calculation 
 with limited sampling can be very untrustworthy . in principal , nontempered metad can also provide converged 
 pmf estimates with two levels of averaging : time averaging over a 
 period of the simulation after it becomes diffusive and replica averaging 
 over the time - averaged bias potential .
however , the diffusive motion 
 itself is difficult to realize due to the limited sampling in ntp 
 hydrolysis . to the best of our knowledge
, this work is the first time 
 that diffusive metad behavior has been demonstrated for a qm / mm simulation 
 of an ntp hydrolysis reaction . the enhanced sampling methods 
 described in this paper , ttmetad 
 and mbmetad , have been implemented in a publicly available fork of 
 the latest version of plumed2 and are free to use .
we are currently 
 working on implementing ttmetad and mbmetad into the public version 
 of the most recent plumed .
both ttmetad and mbmetad are easy to understand 
 and apply , and , compared to a conventional nontempered metad or wtmetad 
 simulation , the extra input parameters that ttmetad and mbmetad require 
 have robust defaults and are straightforward to incorporate .
no specialized 
 or separate mathematical libraries are required by the new features 
 of these methods , thus any researcher already using metad and plumed2 
 can begin using ttmetad and mbmetad simply by reinstalling the new 
 fork of plumed2 with any compatible md engine .
an example mbmetad 
 input file ( used in our simulation ) is provided in the supporting information . to conclude , in 
 this paper we have demonstrated the combined advantages 
 of ttmetad and mbmetad simulation for demanding qm / mm md free energy 
 simulations of the atp hydrolysis process in actin .
in turn , the calculated 
 ( and demonstrably reasonably well converged ) difference in free energy 
 barriers for g - actin versus f - actin is found to be in good agreement 
 with the experimental value .
our work therefore serves as an example 
 for utilizing this approach in simulating complex ntp hydrolysis reactions 
 in proteins , which are some of the most challenging of all biochemical reactions .
additional 
 simulations along these lines will be the target of future research , 
 as well as additional mechanistic analysis .
 the 
 protein mediated hydrolysis of nucleoside triphosphates such 
 as atp or gtp is one of the most important and challenging biochemical 
 reactions in nature . 
 the chemical environment ( water structure , catalytic 
 metal , and amino acid residues ) adjacent to the hydrolysis site contains 
 hundreds of atoms , usually greatly limiting the amount of the free 
 energy sampling that one can achieve from computationally demanding 
 electronic structure calculations such as qm / mm simulations . 
 therefore , 
 the combination of qm / mm molecular dynamics with the recently developed 
 transition - tempered metadynamics ( ttmetad ) , an enhanced sampling method 
 that can provide a high - quality free energy estimate at an early stage 
 in a simulation , is an ideal approach to address the biomolecular 
 nucleoside triphosphate hydrolysis problem . in this work 
 the atp hydrolysis 
 process in monomeric and filamentous actin is studied as an example 
 application of the combined methodology . 
 the performance of ttmetad 
 in these demanding qm / mm simulations is compared with that of the 
 more conventional well - tempered metadynamics ( wtmetad ) . 
 our results 
 show that ttmetad exhibits much better exploration of the hydrolysis 
 reaction free energy surface in two key collective variables ( cvs ) 
 during the early stages of the qm / mm simulation than does wtmetad . 
 
 the ttmetad simulations also reveal that a key third degree of freedom , 
 the o  h bond - breaking and proton transfer from the lytic water , 
 must be biased for ttmetad to converge fully . to perturb the ntp hydrolysis 
 dynamics to the least extent and to properly focus the metad free 
 energy sampling , we also adopt here the recently developed metabasin 
 metadynamics ( mbmetad ) to construct a self - limiting bias potential 
 that only applies to the lytic water after its nucleophilic attack 
 of the phosphate of atp . with these new , state - of - the - art enhanced 
 sampling metadynamics techniques 
 , we present an effective and accurate 
 computational strategy for combining qm / mm molecular dynamics simulation 
 with free energy sampling methodology , including a means to analyze 
 the convergence of the calculations through robust numerical criteria . 
Introduction
Methodology
Results
Conclusions