dynamic 
 biomolecules often undergo large - scale structural changes 
 and visit distinct conformational states during their biological function .
it is of great biological and pharmaceutical 
 interest to characterize their structures and conformational transition 
 pathways .
ideally , detailed free energy landscapes are sought to understand 
 the functional mechanisms of biomolecules in a quantitative manner . 

however , due to large energy barriers , conformational transitions 
 of biomolecules usually take place on time scales of milliseconds 
 or even longer .
this has presented a grand challenge 
 for computational molecular dynamics ( md ) simulations that are limited 
 to typically hundreds - of - nanoseconds to tens - of - microseconds . to address this challenge ,
biasing simulation 
 methods have been 
 found useful in enhanced sampling and free energy calculation of the 
 biomolecules .
these methods include umbrella sampling , conformational flooding , metadynamics , adaptive biasing force ( abf ) calculations , and orthogonal space sampling , etc . during the simulations , a potential or force bias 

is applied along certain reaction coordinates ( or collective variables ) 
 to facilitate the biomolecular conformational transitions across high - energy 
 barriers .
typical reaction coordinates include atom distances , torsional 
 dihedrals , root - mean - square deviation ( rmsd ) relative to a reference 
 configuration , eigenvectors generated from the principal component 
 analysis , and so on .
the definition of 
 the reaction coordinates , however , often requires expert knowledge 
 of the studied systems .
furthermore , the predefined reaction coordinates 
 largely place constraints on the pathway and conformational space 
 to be sampled during the biasing simulations .
it often leads to slow 
 convergence of the simulations when important reaction coordinates 
 are missed during the simulation setup . accelerated molecular dynamics ( amd ) is an enhanced sampling 
 technique 
 that works often by adding a non - negative boost potential to smoothen 
 the biomolecular potential energy surface .
the boost potential , v , decreases the energy barriers and thus accelerates transitions 
 between the different low - energy states . with this
, amd is able to sample distinct biomolecular conformations 
 and rare barrier - crossing events that are not accessible to conventional 
 md ( cmd ) simulations . unlike the previously mentioned biasing simulation 
 methods
, amd does not require predefined reaction coordinate(s ) , which 
 can be advantageous for exploring the biomolecular conformational 
 space without a priori knowledge or restraints .
amd has been successfully 
 applied to a number of biological systems , ( 10 ) and hundreds - of - nanosecond amd simulations have been shown to capture 
 millisecond - time - scale events in both globular and membrane proteins . while amd has been demonstrated to greatly 
 enhance conformational 
 sampling of biomolecules , it suffers from large energetic noise during 
 reweighting .
the boost potential applied 
 in amd simulations is typically on the order of tens - to - hundreds of 
 kilocalories per mole , which is much greater in magnitude and wider 
 in distribution than that of other biasing simulation methods that 
 make use of predefined reaction coordinates ( e.g. , several kilocalories 
 per mole ) .
it has been a long - standing problem to accurately reweight 
 amd simulations and recover the original free energy landscapes , especially 
 for large proteins .
our recent studies 
 showed that when the boost potential follows near - gaussian distribution , 
 cumulant expansion to the second order provides improved reweighting 
 of amd simulations compared with the previously used exponential average 
 and maclaurin series expansion reweighting methods .
the reweighted free energy profiles are in good agreement 
 with the long - time - scale cmd simulations as demonstrated on alanine 
 dipeptide and fast - folding proteins .
however , 
 such improvement is limited to rather small systems ( e.g. , proteins 
 with less than 35 amino acid residues ) . in simulations of larger systems , the boost potential exhibits 
 significantly wider distribution and
, a gaussian accelerated molecular dynamics ( gamd ) approach 
 is presented to reduce the energetic noise for simultaneous unstrained 
 enhanced sampling and free energy calculation of biomolecules , even 
 for large proteins .
gamd makes use of harmonic functions to construct 
 the boost potential that is adaptively added to the biomolecular potential 
 energy surface .
a minimal set of simulation parameters is dynamically 
 adjusted to control the magnitude and distribution width of the boost 
 potential . as such
, the resulting boost potential follows gaussian 
 distribution and allows for accurate reweighting of the simulations 
 using cumulant expansion to the second order . in this study ,
gamd 
 is demonstrated on unconstrained simulations of the alanine dipeptide , 
 chignolin folding , and ligand binding to the t4-lysozyme .
gaussian accelerated 
 molecular dynamics enhances the conformational sampling of biomolecules 
 by adding a harmonic boost potential to smoothen the system potential 
 energy surface ( figure 1 ) .
consider a system with n atoms at positions   = { 1, ... ,n } .
when the 
 system potential v ( ) is lower than 
 a threshold energy , e , a boost potential is added 
 as1where k is 
 the harmonic force 
 constant .
the modified system potential , v * ( ) = v ( ) + v ( ) , is given 
 by2otherwise , when the system potential 
 is above 
 the threshold energy , i.e. , v ( )  e , the boost 
 potential is set to zero and v * ( ) = v ( ) .
when the threshold energy is set to the maximum potential 
 ( e = vmax ) , the system 
 potential energy surface is smoothened by adding a harmonic boost 
 potential that follows gaussian distribution .
the coefficient k0 in the range of 01 determines the 
 magnitude of the applied boost potential . with greater k0 ,
higher boost potential is added to the original energy 
 surface in conventional molecular dynamics ( cmd ) , which provides enhanced 
 sampling of biomolecules across decreased energy barriers . in order to smoothen the potential energy surface 
 for enhanced 
 sampling ,
first , for any two arbitrary potential values v1 ( ) and v2 ( ) found on the 
 original energy surface , 
 if v1 ( ) < v2 ( ) , v should be a 
 monotonic function that does not change the relative order of the 
 biased potential values ; i.e. , v1 * ( ) < v2 * ( ) . by replacing v * ( ) with eq 2 and isolating e , we then obtain3second , 
 if v1 ( ) < v2 ( ) , the potential 
 difference observed on 
 the smoothened energy surface should be smaller than that of the original ; 
 i.e. , v2 * ( )  v1 * ( ) < v2 ( )  v1 ( ) .
similarly , 
 by replacing v * ( ) with eq 2 , we can derive4with vmin  v1 ( ) < v2 ( )  vmax , we need to set the threshold energy e in the following 
 range by combining eqs 3 and 4:5where vmin and vmax are the system 
 minimum and maximum potential 
 energies . to ensure that eq 5 is valid , vmax  vmin + ( 1/k ) and k has to satisfy the following:6we define k  k0(1/(vmax  vmin ) ) ; then 0 < k0  1 . as illustrated in figure 1 , k0 determines the magnitude 
 of the applied boost potential . with greater k0 ,
higher boost potential is added to the potential energy 
 surface , which provides enhanced sampling of biomolecules across decreased 
 energy barriers .
third , the standard deviation of v needs 
 to be small enough ( i.e. , narrow distribution ) to ensure accurate 
 reweighting using cumulant expansion to the second order:7where vav and 
 v are the average and standard 
 deviation of the system potential energies and v is the standard deviation of v with 0 as a user - specified upper limit ( e.g. , 
 10kbt ) for accurate reweighting .
provided eq 5 gives 
 the range of threshold energy e , when e is set to the lower bound e = vmax , we plug in e and k and obtain8we define the right - hand side in eq 8 as k0 = ( 0/v)((vmax  vmin)/(vmax  vav ) ) . for efficient 

enhanced sampling with the highest possible acceleration , k0 can then be set to its upper bound as9the greater v that is obtained from the original potential energy surface 
 ( particularly for large biomolecules ) , the smaller k0 may be applicable to allow for accurate reweighting . 
 alternatively ,
when the threshold energy e is set 
 to its upper bound e = vmin + ( 1/k ) according to eq 5 , we plug in e and k in eq 7 and 
 obtain10we define the 
 right - hand side in eq 10 as k0  
 ( 1  ( 0/v ) ) ( ( vmax  vmin)/(vav  vmin ) ) .
when k0  0 with v  0 , k0 can be theoretically set to any value between zero and 
 1 , although k0 = 1 is applied in this 
 case for the current implementation of gamd ( see appendix a ) . note that a smaller k0 will give higher threshold energy e , but smaller 
 force constant k. when 0 < k0  
 1 , k0 can be set to either k0 for 
 the highest threshold energy e or its upper bound 
 1.0 for the greatest force constant k. in this case k0 = k0 is applied in the current gamd 
 implementation . when k0 >
1 , we have to lower the 
 threshold 
 energy e to ensure that 0 < k0  1 ; e.g. , e = vmax for default k0 = 1 . given e and k0 , we 
 can calculate the boost potential as11similar 
 to amd , gamd provides 
 options to add only the total potential boost vp , only dihedral potential boost vd , or the dual potential boost ( both vp and vd ) .
the dual - boost 
 simulation generally provides higher acceleration than the other two 
 types of simulations for enhanced sampling .
the simulation parameters comprise of the threshold energy values 
 and the effective harmonic force constants , k0p and k0d for the total and dihedral 
 potential boost , respectively . to characterize the extent to 
 which v follows 
 gaussian distribution ,
its distribution anharmonicity ,  , is 
 calculated as was done previously:12where v is 
 dimensionless 
 as divided by kbt with kb and t being the boltzmann 
 constant and system temperature , respectively , and smax = ( 1/2 ) ln(2ev2 ) is the maximum entropy of v .
reweighting 
 by approximating the exponential average term with cumulant expansion 
 to the second order is able to accurately recover the original free 
 energy landscape ( appendix b ) .
as  
 increases , the v distribution becomes less 
 harmonic and the reweighted free energy profile obtained from cumulant 
 expansion to the second order deviates from the original . as demonstrated 
 on the alanine dipeptide , chignolin , and t4-lysozyme in this study , 
 a gaussian distribution of the boost potential is normally achieved 
 in the gamd simulations .
the anharmonicity of v distribution serves as an indicator of the enhanced sampling convergence 
 and accuracy of the reweighted free energy .
gamd is currently implemented 
 in the gpu version of amber 12 ( see appendix a for implementation details ) .
simulations 
 of the alanine dipeptide , chignolin , and t4-lysozyme were performed 
 using the amber ff99sb force field on gpus .
, the ligand 
 benzene was removed from the x - ray crystal structure of the leu99ala 
 mutant ( protein data bank ( pdb ) , 181l ) .
another four benzene molecules were 
 placed in the bulk solvent at least 40  away from the ligand - binding 
 site in the starting configuration . by solvating the structures in 
 a tip3p water box that extends 810 
  from the solute surface , the alanine dipeptide system contained 
 630 water molecules , 2,211 waters for chignolin , and 9,011 waters 
 for t4-lysozyme .
the total number of atoms in the three systems are 
 1,912 ; 6,773 ; and 29,692 for alanine dipeptide , chignolin , and t4-lysozyme , 
 respectively ( table 1 ) .
bonds containing hydrogen atoms were restrained with the 
 shake algorithm , and a 2 fs time step 
 was used .
weak coupling to an external temperature and pressure bath 
 was used to control both temperature and pressure .
the electrostatic interactions were calculated using the 
 pme ( particle mesh ewald summation ) with 
 a cutoff of 8.0  for long - range interactions .
the three 
 systems were initially minimized for 2,000 steps using 
 the conjugate gradient minimization algorithm , and then the solvent 
 was equilibrated for 50 ps in an isothermal  isobaric ( npt ) ensemble with the solute atoms fixed .
another minimization 
 was performed with all atoms free , and the systems were slowly heated 
 to 300 k over 500 ps .
final system equilibration was achieved by a 
 200 ps isothermal  isovolumetric ( nvt ) and 
 400 ps npt run to ensure that the water box of the 
 simulated systems had reached the appropriate density .
in the 
 present study , the system threshold energy is set as e = vmax for all gamd simulations . 

the maximum , minimum , average , and standard deviation values of the 
 system potential ( vmax , vmin , vav , and v ) were obtained from an initial 2 ns nvt simulation with no boost potential . for optimal acceleration , 
 the greatest 0 and k0
were determined through short testing simulations ( e.g. , 2 ns ) : 
 with increasing 0 , either k0 calculated using eq 9 reaches 1.0 ( highest acceleration level ) or v reaches 10kbt ( the upper limit of the v distribution 
 width for accurate reweighting ) .
each gamd simulation proceeds with 
 a 2 ns equilibration run , followed by production simulations .
dihedral potential boost  
 only to identify the optimal 0 and k0 acceleration parameters ( supporting 
 information tables s1s3 ) .
long - time dual - boost gamd 
 simulations were then obtained for analysis , including three independent 
 30 ns simulations of alanine dipeptide and three independent 300 ns 
 simulations of chignolin .
complete binding of benzene 
 to the target ligand - binding site was observed in one of the five 
 simulations . even when the simulation was extended to 1,800 ns ,
time courses of dihedral angles , 
 rmsd , radius of gyration ( rg ) , and residue 
 distances in the amber simulation trajectories were analyzed using 
 the cpptraj tool .
particularly , 
 the backbone dihedral angles  and  were calculated for 
 alanine dipeptide ( figure 2a ) .
for chignolin , the rg and 
 rmsd of simulation frames relative to the pdb native structure ( figure 3a ) were calculated 
 for the protein c atoms with the terminal residues 
 gly1 and gly10 excluded . for ligand binding to the t4-lysozyme ( figure 5a )
, the symmetry - corrected 
 rmsd of benzene was obtained by generating six symmetrically imaged 
 reference benzene configurations from the 181l crystal structure , calculating the rmsds 
 for the diffusing benzene molecules in each frame after aligning the 
 protein c - terminal domain ( residues 80160 ) and then extracting 
 the minimum value of the calculated rmsds .
moreover , the number of 
 protein atoms found within 5  of benzene ( ncontact , only heavy atoms are considered ) was calculated 
 using pbwithin in vmd that accounts for the periodic 
 boundary conditions .
the pyreweighting toolkit was used 
 to reweight the gamd simulations for calculating the pmf profiles 
 and to examine the boost potential distributions .
two - dimensional 
 ( 2d ) pmf profiles were computed for backbone dihedrals (  ,  ) 
 in alanine dipeptide . a bin size of 6 is selected to balance 
 between reducing the anharmonicity and increasing the bin resolution 
 as discussed earlier .
two - dimensional 
 pmf profiles were also constructed using ( rmsd , rg ) for chignolin with a bin size of ( 1.0  , 1.0  ) . 
 for benzene binding to the t4-lysozyme
, 2d pmf was constructed using 
 ( ligand rmsd , ncontact ) with a bin size 
 of ( 1.0  , 5 ) . when the number of simulation frames within a 
 bin is lower than a certain limit ( i.e. , cutoff ) , the bin is not sufficiently 
 sampled and thus is excluded for reweighting .
the cutoff can be determined 
 by iteratively increasing it until the minimum position of the pmf 
 profile does not change .
the final cutoff 
 was set as 10 , 50 , and 1000 for reweighting of gamd simulations on 
 alanine dipeptide , chignolin , and t4-lysozyme , respectively .
gamd provides enhanced sampling of conformational transitions in 
 alanine dipeptide , chignolin folding , and ligand binding to the t4-lysozyme . 

furthermore , the boost potential applied in the present gamd simulations 
 follows gaussian distribution , which allows accurate reweighting using 
 cumulant expansion to the second order and recovery of the original 
 biomolecular free energy landscapes , even for proteins as large as 
 the t4-lysozyme .
notably , hundreds - of - nanoseconds gamd simulations 
 are able to capture complete folding of chignolin and ligand binding 
 of benzene to the t4-lysozyme that take place on significantly longer 
 time scales . in order to balance between achieving 
 the highest acceleration ( large v ) and ensuring 
 accurate reweighting ( small enough standard deviation , v ) , short gamd simulations ( 2 ns ) were 
 tested on alanine dipeptide to search for optimal acceleration parameters . 

for the total potential boost , when 0p was adjusted 
 from 1.0 to 2.0 , the resulting k0p calculated 
 using eq 9 was increased 
 from 0.21 to 1.0 , accompanied by increases of the vp standard deviation from 1.03 to 1.75 kcal / mol and the 
 vp average from 2.36 to 3.85 kcal / mol 
 ( see supporting information table s1a ) . 
 after k0p reached its maximum 1.0 ,
slight 
 changes were observed when 0p was further increased 
 from 2.0 to 3.0 . for the dihedral potential boost ,
the calculated k0d reached 1.0 when 0d was 
 increased to 3.0 ( supporting information table s1b ) .
therefore , ( 0p,0d ) were set to ( 3.0 , 3.0 ) for the production dual - boost gamd simulations 
 of alanine dipeptide , which enables the highest acceleration as well 
 as accurate reweighting with both k0p and k0d equal to 1.0 ( supporting 
 information table s1c ) .
as shown in figure 2b , the boost potential v applied in gamd 
 simulation of the alanine dipeptide follows gaussian distribution . 

the average and standard deviation of v are 10.9 and 2.9 kcal / mol , respectively ( table 1 ) .
figure 2c plots the 2d pmf of backbone dihedrals 
 (  ,  ) that are obtained by reweighting the three 30 ns 
 gamd simulations combined using cumulant expansion to the second order . 

the reweighted pmf is able to recover five energy minimum wells that 
 are centered around ( 162 , 12 ) and ( 72 , 
 12 ) for the right - handed  helix ( r ) , ( 48 , 18 ) for the left - handed  helix 
 ( l ) , ( 156 , 162 ) for the -sheet , 
 and ( 66 , 156 ) for the polyproline ii ( pii ) conformation .
the corresponding minimum free energies are estimated 
 as 0 , 1.40 , 1.07 , 2.94 , and 4.27 kcal / mol , respectively .
this is in 
 good agreement with the pmf results of an exceptionally long 1000 
 ns cmd simulation as presented earlier .
additionally , the distribution anharmonicity of vof frames clustered in each bin of the 2d pmf exhibits values smaller 
 than 0.10 in the low - energy regions ( figure 2d ) , which suggests sufficient sampling for 
 reweighting using cumulant expansion to the second order . greater 
 anharmonicity
is found in high - energy regions ( especially the energy 
 barriers ) , for which low boost potential is applied and less sampling 
 is normally achieved .
therefore , the anharmonicity of v distribution appears to be a good indicator for sufficiency 
 of the enhanced sampling and accuracy of the reweighted free energy .
( a ) scheme 
 representation of backbone dihedrals  and  
 in alanine dipeptide .
( c ) 2d potential of mean 
 force ( pmf ) of backbone dihedrals (  ,  ) calculated from 
 three 30 ns gamd simulations combined using cumulant expansion to 
 the second order .
the low - energy wells are labeled corresponding to 
 the right - handed  helix ( r ) , left - handed 
  helix ( l ) , -sheet (  ) , and 
 polyproline ii ( pii ) conformations .
( d ) distribution anharmonicity 
 of v of frames found in each bin of the pmf 
 profile . for chignolin
, short 2 ns testing 
 gamd simulations showed that k0p reached 
 the 1.0 maximum for the highest total potential boost when 0p was increased to 3.0 ( supporting information table s2a ) .
for the dihedral potential boost , the simulation crashed 
 when 0d was increased to 0.9 and stable simulation 
 was achieved only when 0d was increased to 0.6 ( supporting information table s2b ) .
thus , ( 0p,0d ) were set to ( 3.0 , 0.6 ) for the production 
 dual - boost gamd simulations of chignolin ( supporting 
 information table s2c ) .
the resulting boost potential follows 
 gaussian distribution and the calculated  = 9.22  10 ( figure 3b ) .
the average 
 and standard deviation of v are 9.8 and 2.4 
 kcal / mol , respectively ( table 1 ) .
folding of chignolin captured in gamd simulations : ( a ) comparison 
 of simulation - folded chignolin ( blue ) with the pdb ( 1uao ) native structure 
 ( red ) that exhibits 0.2  rmsd , ( b ) distribution of the boost 
 potential v , ( c ) 2d ( rmsd , rg ) pmf calculated by reweighting the three 300 ns gamd 
 simulations combined , and ( d ) the distribution anharmonicity of v of frames found in each bin of the pmf profile . started from an extended conformation , 
 gamd simulations were able 
 to capture complete folding of chignolin into its native structure 
 ( supporting information movie s1 ) .
the 
 rmsd obtained between the simulation - folded chignolin and nmr experimental 
 native structure ( pdb , 1uao ) reaches a minimum of 0.2  ( figure 3a ) .
using the protein rmsd 
 relative to the pdb native structure and rg , a 2d pmf profile was calculated by reweighting the three 300 ns 
 gamd simulations combined ( figure 3c ) .
the reweighted pmf allows us to identify three 
 distinct low - energy conformational states : the folded (  f  ) , 
 unfolded (  u  ) , and an intermediate (  i  ) . 

the folded state corresponds to the global energy minimum at ( 1.0 
  , 4.0  ) , the unfolded state is 3.68 kcal / mol higher in 
 a local - energy well centered at ( 6.0  , 7.0  ) , and the 
 intermediate in a 3.06 kcal / mol free energy well centered at ( 4.0 
  , 5.5  ) .
the energy barrier for chignolin folding between 
 the unfolded and folded states is 4.0 kcal / mol , and 3.5 
 kcal / mol for transitions between the intermediate and folded states 
 ( figure 3c ) . in comparison , 
 when the three 300 ns gamd simulations were analyzed separately , the 
 reweighted pmf profiles exhibit significant differences ( supporting information figure s2 ) .
whereas the 
 three folded , intermediate , and unfolded low - energy states were captured 
 in sim2 , the unfolded and intermediate states were not sufficiently 
 sampled ( no converged low - energy wells ) in sim1 and sim3 , respectively . 

figure 3d plots the distribution anharmonicity of v for frames found in each bin of the 2d pmf as shown in figure 3c .
the anharmonicity 
 exhibits values smaller than 0.05 in the simulation sampled conformational 
 space , suggesting that v achieves sufficient 
 sampling for reweighting using cumulant expansion to the second order .
compared with the native structure , the unfolded chignolin exhibits 
 an extended conformation without proper secondary structure formed 
 in the protein backbone as shown in figure 4a .
the intermediate conformation observed 
 during folding of chignolin is characterized by hydrophobic interactions 
 between the pro4 and trp9 residue side chains and a turn at thr8 . 

such conformation was also observed in previous microsecond - time - scale 
 cmd simulations . in the simulation - derived 

folded state , the protein residue side chains exhibit closely similar 
 conformations as in the nmr native structure ( figure 4c ) . residues tyr2 and pro4 form hydrophobic 
 interactions between their side chains .
hydrophilic residues asp3 , 
 glu5 , thr6 , and thr8 expose side chains to the solvent .
structures 
 of chignolin observed during gamd simulations corresponding 
 to the ( a ) unfolded (  u  ) , ( b ) intermediate (  i  ) , 
 and ( c ) folded (  f  ) states ( blue ) aligned to the pdb 
 native structure ( red ) .
residues including tyr2 , asp3 , pro4 , glu5 , 
 thr6 , thr8 , and trp9 are represented by sticks .
finally , a 2d free energy profile 
 of chignolin was calculated using 
 the rmsd of the protein c atoms relative to the 
 native pdb structure and the backbone dihedral  of residue 
 gly7 ( supporting information figure s3 ) . 

this allows for identification of a misfolded (  m  ) 
 low - energy state , in addition to the f , i , and u states as observed 
 in figure 3c .
the m 
 conformation of chignolin shows 180 rotation of the 
 c - terminal strand about its long axis relative to the native pdb structure . 

the free energy of the misfolded chignolin is found to be 1.89 kcal / mol 
 greater than that of the folded state .
this is consistent with previous 
 findings that the misfolded chignolin was observed with higher probability 
 in simulations with the amber ff99sb force field , which is used in the present study .
short testing gamd simulations 
 of the t4-lysozyme showed that k0p reached 
 the 1.0 maximum for the highest total potential boost when 0p was increased to 3.0 ( supporting information table s3a ) .
for the dihedral potential boost , the simulation crashed 
 when 0d was increased to 5.0 ( supporting information table s3b ) .
stable simulation was achieved 
 when 0d was increased to 4.0 and the calculated k0d was 0.35 . with further testing for stable 
 dual - boost gamd simulation , ( 0p,0d ) were set to ( 3.0 , 4.0 ) for the final production simulations of 
 the t4-lysozyme ( supporting information table s3c ) .
gamd captured complete binding of benzene to the 
 deeply buried ligand - binding cavity in the leu99ala t4-lysozyme within 
 100 ns in one of the five independent 800 ns simulations ( supporting information movie s2 ) .
benzene remained 
 bound in the ligand - binding site even when the simulation was extended 
 to 1,800 ns . by aligning the c - terminal domain ( residues 80160 ) 
 of the t4-lysozyme , the rmsd of the diffusing benzene molecules relative 
 to the bound pose in the 181l x - ray crystal structure reaches a minimum of 0.1  
 ( figure 5a ) .
the boost potential applied during the 1,800 ns 
 gamd simulation follows gaussian distribution , and its  = 1.39 
  10 ( figure 5b ) .
the average and standard deviation of v are 36.5 and 4.7 kcal / mol , respectively ( table 1 ) .
although the v average values exhibit variations between five independent 
 simulations , the v standard deviations are 
 closely similar to each other provided that ( 0p,0d ) were set to ( 3.0 , 4.0 ) ( supporting 
 information table s3c ) . using the rmsd of benzene relative 
 to the bound pose and the number of protein heavy atoms that are within 
 5  of benzene ( ncontact ) ,
a 2d pmf 
 profile was calculated by reweighting the 1,800 ns gamd simulation 
 ( figure 5c ) .
the reweighted pmf allows us to identify three 
 distinct low - energy states : the unbound (  u  ) , intermediate 
 (  i  ) , and bound (  b  ) states .
the bound 
 state corresponds to the global energy minimum located at (0 
  , 30 ) , the unbound state in a local - energy well centered at 
 (33  , 0 ) , and the intermediate centered at (11 
  , 20 ) .
it is important to note that since the complete binding 
 of benzene to the target ligand - binding site was observed only once , 
 the calculated binding free energy between the bound and unbound states 
 is subject to the error of limited sampling .
nevertheless , benzene 
 visits the intermediate site many times during the 1800 ns gamd simulation 
 with the ligand rmsd decreased to 11  ( supporting information figure s2a ) .
repeated sampling of the 
 intermediate state was observed in the other four 800 ns gamd simulations 
 as well ( supporting information figure 
 s2 ) , for which a local - energy well appears around ( 11.0  , 20 ) 
 in the 2d pmf profiles ( supporting information figure s3 ) .
the relative free energy between the intermediate and 
 unbound states is estimated to be 0.53  0.46 kcal / mol from pmf 
 profiles of the five gamd simulations .
furthermore , benzene was observed 
 to bind another intermediate 2 (  i2  ) site that is located 
 in the pocket formed by the hinge c helix and the 
 b helix from the n - terminal domain ( supporting information figure s2 ) .
a corresponding 
 local energy well of the i2 state appears in the calculated 2d pmf 
 profiles .
figure 5d 
 plots the v distribution anharmonicity ,  , 
 for frames found in each bin of the 2d pmf .
it exhibits relatively 
 large values in the high - energy regions ( less sampling ) , notably the 
 energy barrier between the intermediate and bound states . the ligand 
 entry from the intermediate to the bound state
is thus suggested to 
 be the rate - limiting step for benzene binding . in comparison ,  
 exhibits values smaller than 0.01 in the energy well regions , suggesting 
 that v achieves sufficient sampling for reweighting 
 using cumulant expansion to the second order .
binding of benzene to 
 the leu99ala t4-lysozyme simulated via gamd : 
 ( a ) comparison of simulation - derived complex structure that captures 
 benzene binding ( blue ) with 0.1  ligand rmsd relative to the 
 181l pdb structure ( red ) , ( b ) distribution of the boost potential 
 v , ( c ) 2d ( ligand rmsd , ncontact ) pmf calculated by reweighting the 1,800 ns gamd 
 simulation , and ( d ) distribution anharmonicity of v of frames found in each bin of the free energy profile .
a complete binding pathway of benzene that was 
 observed in the 
 gamd simulation is shown in figure 6a .
benzene diffuses from the bulk solvent to the protein 
 surface formed by the d and g helices 
 and then to the target ligand - binding site in the protein c - terminal 
 domain .
figure 6b depicts 
 the unbound pose of a benzene molecule that is located far away from 
 the ligand - binding site in the 181l x - ray crystal structure . in the intermediate 
 state ,
benzene interacts with residues lys83 , pro86 , and val87 from 
 the d helix and the thr115 , thr119 , and gln122 residues 
 from the g helix ( figure 6c ) . in the bound pose , benzene is superimposable 
 with the ligand cocrystallized in the 181l crystal structure .
it forms hydrophobic 
 interactions with residues ile78 , leu84 , tyr88 , val87 , leu91 , val111 , 
 leu118 , and leu121 in the deeply buried protein cavity ( figure 6d ) .
supporting 
 information figure s5 shows a transient snapshot observed during 
 benzene binding between the intermediate and bound poses .
benzene 
 appears to slide into the ligand - binding cavity , interacting with 
 residues gln81 , lys83 , leu84 , val111 , phe114 , thr115 , and leu118 from 
 the d , f , and g helices . ( a ) pathway of benzene binding to the t4-lysozyme observed during 
 the gamd simulation .
( b  d ) the ( b ) unbound (  u  ) , 
 ( c ) intermediate (  i  ) , and ( d ) bound (  b  ) 
 poses of the protein  ligand complex ( blue ) with the protein 
 c - terminal domain ( residues 80160 ) aligned to the pdb native 
 structure ( red ) .
the protein and benzene are represented by ribbons 
 and spheres , respectively , and they are colored by blue for the simulation 
 structure while red for the pdb native structure , except that in a 
 the simulated benzene is represented by lines and colored by simulation 
 time in a bwr color scale .
residues with heavy atoms found within 
 3  of benzene are represented by sticks .
by adaptively adding a harmonic boost potential 
 to smoothen the 
 system energy surface , gamd provides both unconstrained enhanced sampling 
 and free energy calculation of biomolecules .
important statistical 
 properties of the system potential , such as the average , maximum , 
 minimum , and standard deviation values , are used to calculate the 
 simulation acceleration parameters , particularly the threshold energy , e , and force constant , k0 . for alanine dipeptide
, both k0p and k0d can be increased to the 1.0 maximum for the 
 greatest possible boost of the total and dihedral potential energies . 

the resulting standard deviation of v that 
 follows gaussian distribution is 2.9 kcal / mol , which allows for accurate 
 reweighting using cumulant expansion to the second order .
notably , 
 high - energy regions of the gamd sampled free energy surface ( especially 
 the energy barriers ) exhibit increased anharmonicity compared with 
 the low - energy wells ( figure 2d ) .
thus , the free energy barriers appear unconverged and 
 still suffer from insufficient sampling during the three short 30 
 ns gamd simulations . nevertheless , the gamd reweighted pmf profile 
 is able to recover low - energy states of the system that are in good 
 agreement with those observed in a 1000 ns cmd simulation .
furthermore , gamd is able to fold chignolin 
 from an extended conformation 
 to the nmr native structure in three independent 300 ns simulations . 

while k0p can be increased to 1.0 for 
 the greatest total potential boost , k0d can reach only 0.150.17 for stable simulation .
nevertheless , 
 this seems to provide sufficient sampling of the chignolin - folding 
 process . among the three independent simulations , although the folded , 
 intermediate , and unfolded low - energy states are captured in sim2 , 
 the unfolded and intermediate states are not sufficiently sampled 
 in sim1 and sim3 , respectively ( supporting information figure s2 ) .
all three folded , intermediate , and unfolded 
 states are identified from the resulting reweighted free energy profiles . 
 particularly , the intermediate conformation was also observed in earlier 
 microsecond - time - scale cmd simulations .
gamd appears to achieve better convergence than the previous amd 
 simulations that were not able to distinguish between the intermediate 
 and unfolded states within the same simulation time . finally , gamd captured complete binding of benzene 
 to the ligand - binding 
 site of the t4-lysozyme . in one of the five independent simulations , 
 benzene diffused into the deeply buried ligand - binding cavity within 
 100 ns and remained bound even when
similar to chignolin , while k0p can be increased to the 1.0 maximum for the greatest total potential 
 boost in the t4-lysozyme simulations , k0d can reach only 0.330.35 for the greatest dihedral boost . 

the resulting v standard deviation is 4.7 
 kcal / mol in the final dual - boost gamd simulations .
such narrow distribution 
 of v ensures accurate reweighting using cumulant 
 expansion to the second order .
distinct low - energy unbound , intermediate , 
 and bound states were identified from the reweighted pmf profiles . 

the atomistic gamd simulation also elucidates a highly detailed binding 
 pathway of benzene that diffuses from the bulk solvent to an intermediate 
 site located on the protein surface formed by the d and g helices and then slides into the target 
 ligand - binding cavity through a channel formed by the d , f , and g helices .
the 
 free energy difference between the intermediate and unbound states 
 was found to be small at 0.53  0.46 kcal / mol as estimated from 
 the five independent gamd simulations .
benzene repeatedly visits the 
 intermediate site on the protein surface . in comparison , the ligand 
 entry from protein surface to the deeply buried protein cavity appears 
 to be the rate - limiting step for complete benzene binding .
it is important 
 to note that the complete ligand binding was not observed in the four 
 800 ns gamd simulations , suggesting that the present gamd simulations 
 still suffer from insufficient sampling of the ligand entry process 
 and the reweighted free energy profiles remain unconverged ( figure 5c and supporting information figure s5 ) .
this is also 
 indicated by the increased anharmonicity corresponding to the free 
 energy barrier between the intermediate and bound states as shown 
 in figure 5d . nevertheless , 
 our gamd simulation captured a binding pathway of benzene to the t4-lysozyme . 

the ligand entry site is indeed adjacent to the mobile f helix ( residues 108113 ) , which has been suggested 
 earlier based on the finding that the 
 f helix exhibits increased b - factors in the leu99ala 
 complex structures compared to the apo structures . here
, gamd is compared with the original amd , particularly 
 on the 
 performance of smoothening the potential energy surface and energetic 
 reweighting . supporting information table 

s4 summarizes the statistical properties of the boost potential and 
 the original and modified potential energies as obtained from amd 
 and gamd simulations of the alanine dipeptide , chignolin , and t4-lysozyme , 
 with reference to the cmd simulations .
specifically , dual - boost amd 
 simulations of the alanine dipeptide and chignolin obtained from our 
 previous study were used for comparison . 
 for the t4-lysozyme , restarting from the initial configuration of 
 gamd simulations , a 200 ns dual - boost amd simulation was performed 
 using the following acceleration parameters : edihed = vdihed_av + 4nres , dihed = 4nres/5 , etotal = vtotal_av + 0.2natoms , and 
 total = 0.2natoms , where vdihed_av and vtotal_av are the average dihedral and total potential energies calculated 
 from a short 10 ns cmd simulation . as shown in supporting information table s4 , although higher average boost 
 potentials are applied in gamd simulations of the three systems than 
 in the amd simulations , the boost potentials exhibit smaller standard 
 deviations ( i.e. , narrower distribution ) in the gamd simulations , 
 except for the alanine dipeptide .
furthermore , anharmonicity of the 
 boost potential distribution is significantly reduced in the gamd 
 simulations .
relative to cmd , both amd and gamd simulations mostly 
 exhibit smaller standard deviations in the modified dihedral and total 
 potential energies .
notably for the t4-lysozyme , standard deviations 
 of the modified potential energies are significantly smaller in the 
 gamd simulations than in amd .
thus , the potential energy surfaces 
 appear more smoothened in gamd enhanced sampling of the t4-lysozyme . 
 provided the narrower distribution and lower anharmonicity of the 
 boost potential ( supporting information table s4 ) , gamd allows for more accurate approximation of exponential 
 average reweighting factor using cumulant expansion to the second 
 order and thus improved free energy calculation , especially for the 
 t4-lysozyme . in the present simulations of ligand binding to 
 the t4-lysozyme , 

complete binding of benzene was observed only once , largely because 
 the target ligand - binding site is deeply buried in the protein . however , 
 benzene binding to the intermediate site on the protein surface formed 
 by the d and g helices was captured 
 many times , which provides better statistics for free energy calculation . 

it has been suggested that a pmf - based approach is appropriate to 
 calculate binding free energy for ligands ( especially the charged ) 
 that bind to the protein surface .
systems 
 of this type for future applications of gamd might include binding 
 of benzamidine to trypsin and allosteric 
 modulators to the protein surface of the g - protein coupled receptors . in comparison with many enhanced sampling 
 methods such as umbrella 
 sampling , conformational flooding , metadynamics , abf 
 calculations , and orthogonal space sampling , gamd has the advantage of no need to set predefined 
 reaction coordinates .
metadynamics , in particular , is another potential 
 biasing technique that has been widely used to map the free energy 
 landscapes of biomolecules such as protein conformational changes and protein  ligand binding . by monitoring the energy surface of biomolecules during the simulation , 

metadynamics keeps adding small gaussians of potential energies to 
 the low - energy regions .
this will eventually fill the low - energy wells 
 and achieve uniform sampling of the free energy surface along selected 
 reaction coordinates .
the usage of predefined coordinates greatly 
 reduces the complexity of biomolecular simulation problems and facilitates 
 the free energy calculations ( e.g. , significantly lower energetic 
 noise compared with amd simulations ) .
however , it is key to select 
 proper reaction coordinates , which often requires expert knowledge 
 of the studied systems .
construction of biomolecular reaction coordinates 
 or collective variables has thus been one of the main objectives in 
 metadynamics studies .
when important 
 reaction coordinates are missed during the simulation setup , metadynamics 
 simulations may suffer from slow convergence problems as discussed 
 earlier .
furthermore , the predefined reaction coordinates tend to 
 place constraints on the sampled space and pathways .
it seems difficult 
 to identify certain intermediate states during the protein - folding 
 and ligand - binding pathways , e.g. , intermediate 2 as observed during 
 the binding of benzene to the t4-lysozyme . in comparison ,
amd simulations 
 are not constrained by reaction coordinates , but this also leads to 
 much higher energetic noise and presents grand challenge for accurate 
 reweighting to recover the original free energy landscapes of biomolecules . although cumulant expansion to the second order 
 was shown to improve amd reweighting when the boost potential follows 
 near gaussian distribution , such improved 
 reweighting is still limited to small systems such as protein with 
 35 residues . here , by constructing 
 boost potential using a harmonic function that follows gaussian distribution , 
 gamd enables rigorous energetic reweighting through cumulant expansion 
 to the second order , even for simulations of larger proteins ( e.g. , 
 t4-lysozyme ) . with this
first , while 
 the present gamd simulations seem to provide sufficient sampling of 
 the low - energy regions , they appear to remain unconverged in sampling 
 of the high - energy barriers .
this is particularly true for the ligand 
 entry step in the gamd simulation of benzene binding to the t4-lysozyme . 

it is worth recalling that the threshold energy for adding the boost 
 potential is set to its lower bound in the present gamd simulations . 

it is subject to future investigation into whether using the upper 
 bound of the threshold energy will facilitate sampling of the high - energy 
 barriers in gamd simulations .
second , based on a potential biasing 
 approach , gamd mainly accelerates transitions across enthalpic energy 
 barriers .
improvement for its application to systems with high entropic 
 barriers is still in need . on this
regard , gamd can be potentially 
 combined with the parallel tempering and 
 replica exchange algorithms for further 
 enhanced sampling .
particularly , the combination of parallel tempering 
 and metadynamics ( pt - metad ) has been 
 shown to facilitate enhanced sampling of biomolecules over entropic 
 barriers . in summary , without the need to set predefined reaction 
 coordinates , 
 gamd is generally applicable to a wide range of biomolecular systems 
 as demonstrated on protein folding and ligand binding in this study . 
 for systems of increasing size , the upper limit of the v standard deviation
, 0 can be adjusted 
 dynamically to ensure that the distribution width of the applied boost 
 potential is narrow enough for accurate energetic reweighting using 
 cumulant expansion to the second order .
therefore , gamd provides both 
 unconstrained enhanced sampling and free energy calculation for biomolecular 
 simulations .
 a gaussian accelerated molecular 
 dynamics ( gamd ) approach for simultaneous 
 enhanced sampling and free energy calculation of biomolecules is presented . 
 by constructing a boost potential that follows gaussian distribution , 
 accurate reweighting of the gamd simulations 
 is achieved using cumulant 
 expansion to the second order . here , gamd is demonstrated on three 
 biomolecular model systems : alanine dipeptide , chignolin folding , 
 and ligand binding to the t4-lysozyme . without the need to set predefined 
 reaction coordinates , 
 gamd enables unconstrained enhanced sampling 
 of these biomolecules . 
 furthermore , the free energy profiles obtained 
 from reweighting of the gamd simulations allow us to identify distinct 
 low - energy states of the biomolecules and characterize the protein - folding 
 and ligand - binding pathways quantitatively . 
Introduction
Methods
Results
Discussion