Free Energy Calculations using a Swarm-Enhanced Sampling Molecular Dynamics Approach

Free energy simulations are an established computational tool in modelling chemical change in the condensed phase. However, sampling of kinetically distinct substates remains a challenge to these approaches. As a route to addressing this, we link the methods of thermodynamic integration (TI) and swarm-enhanced sampling molecular dynamics (sesMD), where simulation replicas interact cooperatively to aid transitions over energy barriers. We illustrate the approach by using alchemical alkane transformations in solution, comparing them with the multiple independent trajectory TI (IT-TI) method. Free energy changes for transitions computed by using IT-TI grew increasingly inaccurate as the intramolecular barrier was heightened. By contrast, swarm-enhanced sampling TI (sesTI) calculations showed clear improvements in sampling efficiency, leading to more accurate computed free energy differences, even in the case of the highest barrier height. The sesTI approach, therefore, has potential in addressing chemical change in systems where conformations exist in slow exchange.


Introduction
Free energy simulation methods such as free energy perturbation (FEP) and thermodynamic integration (TI) are increasingly important tools in the pursuit of molecular design. [1,2] Despite recent advances, however,the computation of free energies remains challenging for many systems, where the free energy change may depend on multiple distinct substates separated by non-negligible energy barriers. [3] Suitably sampling kinetically distinct but thermodynamically relevant substates is challenging,p articularly where prior knowledge of their existence is lacking. [4] Considerable effort has been invested in tackling the problemo fp seudo-ergodicity in simulation-based free energy calculations. [3,[5][6][7] Ar ange of sampling methodologies have been proposed;f or example, approaches based on TI include the independenttrajectory TI (IT-TI) method, [8] enhancing sampling by coupling to accelerated molecular dynamics (aMD) [9] or replica-exchange-based methods such as RETI. [10] Recently,w ehave explored am ulti-copym olecular dynamics methodd esignedt oi mprove conformational exploration of rugged energy landscapes. [11] Ad evelopment out of the SWARM-MD approach of Huber and van Gunsteren, [12,13] we term this method swarm-enhanced sampling molecular dynamics (sesMD). [11] The sesMDmethod links multiple simulation replicas into as warm, using attractive and repulsive pair potentials acting on dihedrala ngles to promote barrierc rossing into alternative energy minima. Application of sesMD has led to enhanced sampling of the conformationso fs mall-molecule systemsa nd ap rotein kinase, [11] as swarm replicas cooperatively sample ag reater volume of phase space by steering each other across potential energy barriers.
The possibility of harnessing the conformational exploration afforded by swarm-coupled trajectories to improve the accuracy of free energy calculations has been recognised previously. [14] Here, we assess the predictiveq uality of free energy estimates by combining the sesMD approach with adualtopology TI framework. We illustrate the approach of this method, denoted hereafter as sesTI, for the diagnostic case of the butaneto-butane alchemical transition in water.A lkane transformations of this type can suffer from errors in the calculated free energy change, owing to inadequate sampling of the hydrocarbon's internal rotations. [9,[15][16][17] For the transition, we study three butane potentials of growing energy barrierh eight between rotamers to represent increasingly distinct energetically low-lyings ubstates. The sesTI scheme is compared with the TI and IT-TI methods, directly examining the effect of swarm-coupling replicas on their sampling of kineticallys eparated substates.

Theory
According to the TI approach, the Helmholtz free energy change DA of at ransition along ac oordinate l is obtained by using Equation (1): Free energy simulations are an established computational tool in modelling chemical change in the condensed phase. However,s ampling of kinetically distinct substates remains ac hallenge to these approaches. As ar oute to addressing this, we link the methods of thermodynamic integration( TI) and swarm-enhanced sampling molecular dynamics (sesMD), where simulation replicas interact cooperativelyt oa id transitions over energy barriers. We illustrate the approach by using alchemicala lkane transformations in solution, comparing them with the multiple independent trajectory TI (IT-TI) method. Free energy changes for transitions computed by using IT-TI grew increasingly inaccurate as the intramolecular barrier was heightened. By contrast, swarm-enhanced sampling TI (sesTI) calculations showedc lear improvements in sampling efficiency, leadingt om ore accurate computed free energyd ifferences, even in the case of the highest barrierh eight. The sesTI approach,t herefore, has potential in addressing chemical change in systemswhere conformations exist in slow exchange.
where l = 0a nd l = 1r epresent the initial and final states of the transition, respectively,a nd V MM is the molecular mechanical potential energy of the system. Thus, the total differencei n free energy DA can be obtained by using an appropriate quadrature scheme to integrate over l from l = 0t ol = 1t he ensemble averaged @V MM ðlÞ @l values. For the sesTI method, we now consider the computationo fDA accordingt ot he TI approach, but within as esMD framework, that is, for as warm of M simulation replicas a.W ef irstly define the total potential, V tot ,a cting on this swarm within as esMD simulation, as shown in Equation (2): where vector r M is the 3NM dimensionv ector describing the coordinates of N atoms in M replicas. Here, V MM (r M )i st he sum of force field potentials, V MM a ðr a Þ.T he swarm-enhanced sampling (ses) potential, V ses ,isd efined by Equation (3): where A-D are suitably calibrated parameters for attractive (A, B) and repulsive (C, D) terms and d ab rms ðf a ; f b Þ is the rootmean-square dihedrala ngle distance of K dihedrals j between swarm members a and b,namely ðK À1 P K j ðf a j À f b j Þ 2 Þ 1=2 . If we consider the contribution of replica a, V ses a þ V MM a ,a n expectation value of dA/dl for the mutationo fasingleb utane molecule can be obtained from the sesMD ensemble average of @V MM ðlÞ @l at l by applying the approachofT orrie and Valleau [18] to recover aB oltzmann weighted average according to Equation (4): Here, we assume that the swarm replicas are weakly coupled, such that their contributionsa re re-weighted according to V ses a .Afurthera pproximation is possible if one assumes that, given sufficient sampling, the time average over ag iven individual replica will converge to the replica average.A pplying this assumption to the context of our sesMD system leads to the expression in Equation (5): In this work, we compare the performance of sesTI with TI and with the non-interacting multi-replica IT-TI approach. [8] Computational Details MD and TI simulations were performed by using the sander module of the AMBER12 software package; [19] sesMD and sesTI calculations were conducted by using am odified version of the sander module of AMBER11. [20] Three model sets of butane parameters were used, each with zero partial charges on the butane atoms. In the first, b1,t he bonded and van der Waals parameters were taken from the GAFF force field. [21] In systems b2 and b3,t he torsional barrier heights were increased from their GAFF values (Table S1, Figure S1);b utane was solvated in ar ectangular box of 584 TIP3P water molecules. [22] Periodic boundary conditions were employed, with ap article-mesh Ewald (PME) treatment [23,24] and a9c ut-off for short-range non-bonded interactions. All MD calculations used a1fs time step and the SHAKE algorithm [25] was applied to constrain solvent bonds. The temperature and pressure were controlled by using Langevin dynamics, [26] with ac ollision frequency of 2ps À1 and aB erendsen barostat [27] with ac oupling constant of 2ps.
Free energy calculations for the butane-to-butane transition for b1, b2 and b3 employed the dual topology approach [28] and soft-core potentials. [29] For TI, IT-TI and sesTI, as traightforward linear scaling [30] of the mutating groups was used, employing a l path of 13 points (l = 0.01, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 0.95 and 0.99). Each window was equilibrated for 800 ps by using rounds of NPT and NVT dynamics. Subsequently,t he final geometry from this trajectory was replicated 12 times, with initial velocities assigned from a2 98 KM axwell-Boltzmann distribution and equilibrated for 100 ps. For IT-TI calculations, this was followed by 5nsN VT production dynamics at 298 K. For sesTI calculations, as es potential was applied to the central f 0 and f 1 angles of butane ( Figure 1). To explore the differing energy landscapes of systems b1, b2 and b3,w ea pplied distinct sets of (A, B) and (C, D) parameters ( Table 1). The ses potentials, with ar epulsive and longer-range attractive profile were empirically fitted to ensure Figure 1. Butane-to-butane alchemicaltransformation through the dualtopology approach. Active groups (in red) and inactive groups (in green)c oexist and alter as mixing parameter l evolves from state 0t o1 .C orresponding centralC-C-C-C dihedrals f 0, and f 1 for states 0a nd 1a re also shown. www.chemphyschem.org as atisfactory frequency of crossing over the gauche-trans energy barrier for each system. Reflecting the greater barrier height from b1 through to b3,t he magnitude of pre-exponential parameters A and Cc orrespondingly increased (Table 1). Prior to a5ns production simulation of the 12-replica swarm of each l,afurther 200 ps equilibration was applied, where the ses potential was gradually increased from zero to its full strength (Table 1). Coordinates were saved for analysis every 1.0 ps and energies every 0.1 ps.
Free energy estimates were computed by using the approach by Steinbrecher et al. [30,31] to obtain independent samples from the windowing trajectories, based on the autocorrelation time t of dV/ dl.T he standard error of the mean, s SEM ,w as computed as s dV=dl ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi t sim =2t p ,w here s dV/dl is the standard deviation and t sim is the total length of the simulation. For IT-TI, free energy estimates were obtained by considering an average over the free energies computed from each of the 12 replicate simulations. For sesTI, free energies were computed from Equations (4) or (5) above. For the latter,b ootstrap sampling [32,33] from the combined 12 trajectories was applied, using uniform sampling and the Mersenne Twister pseudo-random-number generator mt19937. [34] 1.2 10 5 configurations were obtained for each l.G eometrical analysis of trajectories used cpptraj [35] from the AMBER suite.

Results
We consider the ability of sesTI to accurately calculate the free energy changef or the transformation of n-butanet on-butane in explicit aqueous solvent( Figure 1). The path for this test case alchemical transformation involves mutation of the methyl group on carbon C 1 of butane to ah ydrogen atom, and vice versa for C 3 ( Figure 1). However,r egardless of the path l and force fielde mployed, the free energy difference between these states should be zero. We employ here ad ual topologya pproach to TI;t hus, there is the coexistence of mutating groups in the initial reference butanes tate l = 0w ith potential V 0 ,a nd final target butane state l = 1w ith potential V 1 ,a lbeit withoutinteracting during the MD simulation. The remaining atoms, common to both states, evolve according to the potential ð1 À lÞV 0 þ lV 1 .D ue attention is given to the alchemicale nd points by using soft-core van der Waals potentials forb utane and the divisiono fl into 13 suitablys paced windows.
Nevertheless,t he accuracyo ft he TI calculation also relies on the adequate sampling of conformations relevant to the reference and target states across the whole l-dependentp ath, in particular around C-C-C-C torsions f 0 for l = 0a nd f 1 for l = 1 ( Figure 1). These torsion angles can potentially adopt trans (t)a nd gauche (g + and g À )r otamers. In this study,w ec onsider three different butane potentials, b1, b2 and b3.F or all three systems, the t well lies 1kcal mol À1 lower in energyt han the g + and g À minima( Figure S1). However,t he total potential energy barriers from the t to g + /g À wells (i.e. including torsional and non-bonded terms) are progressively enlarged from approximately 4kcal mol À1 (b1)t o6kcal mol À1 (b2)a nd 8kcal mol À1 (b3), as shown by their rotational profiles( Figure S1). To accentuate the issue of sampling, we initially assign t and g À conformations to butane torsions f 0 and f 1 ,r espectively.C onsequently,t hese systems present increasingly challenging access to thermodynamically relevant states for the butane-tobutanem utation.F or the transformationso ft hese three systems, we comparet he performanceo fs esTI with 1) TI calculations based on DA estimates from the 12 individual MD trajectories and 2) DA from averaging over these 12 TI simulations, namely the IT-TI approach. [8] As an initial indication of the ability of unbiasedM Da nd sesMD to surmount energy barriers in b1, b2 and b3,w ec onsider the 5nso fs imulation at l = 0.01, that is, the l point closest to the reference state. Superposition of equal-spaced snapshots from the 12 unbiased MD trajectories of IT-TI at l = 0.01 indicates that, for b1,a ll three f 0 rotamers of butanea re explored ( Figure 2a); for b2,o nly t and g À (Figure 2b)a re explored and only the initial conformation t for b3 (Figure 2c). By comparison, the l = 0.01 window simulation, using 12-replica sesMD samples shows all three rotamers for all three butanem odels ( Figure2a-c). The three models apply ses potentials of increasing strength from b1 to b3 (Table 1);i ta ppears that the broadest coverage of f 0 space is found for b3 (Figure 2c).
From this preliminarya ssessment of underlying sampling, we now turn to consider the estimates of the free energy Table 1. Set of ses potential parameters used in sesTI calculations for systems b1, b2 and b3.

Model
A www.chemphyschem.org change for the butane-to-butane transformation furnished by TI, IT-TI and sesTI approaches. We label the 12 independent TI calculations as TI-01 to TI-12 (Table 2). For the butane system with the lowest energy barriers, b1,wefind that all 12 TI calculations provide ap redicted DA close to zero, with ar ange of À0.06 to 0.03 kcal mol À1 (Table 2). Standard errors are uniformly 0.01 kcal mol À1 in value. Correspondingly,t he combined IT-TI estimate from these individual TI calculations is À0.01 AE 0.01 kcal mol À1 (Table 2), again in good agreement with the theoreticalv alue of zero. For the application of sesTI in enhancing the backbonet orsions of butane during the alchemical change, the estimated DA is 0.01 AE 0.01 kcal mol À1 when reweighting accordingt oE quation (4), and 0.02 AE 0.01 kcal mol À1 when using Equation (5). Thus, the estimates of both schemes are in agreement with each other and close to zero. Reflecting the rotamer sampling at l = 0.01, considered earlier,t he underlying sampling of f 0 and f 1 for both IT-TI and sesTI simulations appeart oe xplore t, g + and g À minima for b1 across all l (first two columns of Figure 3a)a nd with each replica (first two columns of Figure 3b), although to ag reater extent for sesTI. It also appears that the transitions between wells are frequent in the 12 TI simulations. This is exemplified by exploration of the swarm in the l = 0.05 window (Figure 4a) and sampling across the l range ( Figure 5a); for IT-TI, an average transition frequency of 2.8 ns À1 in f 0 are found for the l window simulations (Table 3). For the swarm-coupled sesTI simulations, this frequency increases five-fold to 15.5 ns À1 (Table 3a nd Figure 4a). This differencei ns ampling frequency between IT-TI and sesTI reflects the rate at which DA estimates converge as af unctiono fs imulation time. The individual TI estimateso fDA converge to comparable values that are within 0.2 kcal mol À1 of each other,a ta round2 .5 ns of MD sampling for each l window (Figure 6a). Thisr educes further to within 0.15 kcal mol À1 at 5ns. Similarly,t he IT-TI average based on these individual DA estimates converges within 3nst o À0.01 kcal mol À1 (Figure 7a); the standard error also appears low and stable at 0.01 kcal mol À1 from 2ns( Ta ble 4). Interest-ingly,b oth sesTI estimates provide average DA values of close to zero when using l window simulation lengths of only 100 ps ( Figure 7a); specifically, DA is 0.06 AE 0.06 kcal mol À1 and 0.01 AE 0.06 kcal mol À1 for Equation (4) and Equation (5), respectively (Table 4).
In the second butanes ystem, b2,t he barrier between t and g + is approximately 2kcal mol À1 higher than for b1.T he 12 individual TI calculations for b2 provide am ore variable prediction of DA,r anging from 0.11t o0 .40 kcal mol À1 (Table 2). Indeed, the latter value is provided by calculation TI-04, which diverges somewhatf rom other replicas after 2nso fM D ( Figure 6b). Nevertheless, it is evident that all replicasa re only gradually approaching the correct DA value and require longer than 5nso fM Dp er window ( Figure 6b). Reflecting this, the overall IT-TI estimate of DA is 0.21 AE 0.02 kcal mol À1 (Table 2). For sesTI,t he computed DA is closer to zero, calculated as À0.03 AE 0.03 kcal mol À1 and 0.05 AE 0.03 kcal mol À1 by using Equations (4) and (5), respectively ( Table 2). The sesTI estimates of DA converge more slowly for b2 compared to b1,b ut appear to stabilisea fter 1.6 ns of simulation for each l (Figure 7b). The standard errors in both re-weighting schemes converge by 2ns( Ta ble 4). The augmented rotational barriers in b2 lead to decreased sampling between rotamers, that is, transition frequencies drop by overaf actor of ten for IT-TI simulations and by at hirdf or sesTI simulations (Table 3a sw ell as Figures 4b and 5b). Consequently,w ithin the 5nsw indow, where 12 TI replicas observe only one change of rotamer in b2 on average,t he swarm of sesTI samples 24 such transitions. Clearly,this restricts the overall coverage of the three rotameric states by IT-TI (Figure 3b and5b) and their relative contributions to DA.T his contrasts with the sesTI window,w here each simulation samples all three rotamers (Figure 3b).
Finally,w ec onsider the b3 model of butane, with af urther 2kcal mol À1 increasei nt!g À barrierh eight. Interestingly,t he DA estimatesf rom the 12 TI simulations agree closely with one another,w ith an arrow range of 0.67-0.72 kcal mol À1 ( Table 2). The corresponding IT-TI estimate is 0.70 AE 0.004 kcal mol À1 .F urthermore, these individual (Figure 6c)a nd average TI estimates of DA (Figure 7c)a nd their standard errors (Table 4) converge very rapidly with the l window simulationl ength. By contrast, the sesTI estimates remain reasonable approximationst oz ero, with values of À0.03 AE 0.09 and 0.11 AE 0.23 kcal mol À1 when using Equations (4) and (5), respectively ( Table 2). As expected, the convergenceo ft hese sesTI estimates of DA anda ssociated standard deviations for b3 is slower than for b1 and b2 (Table 4). Thus, the DA computed by using Equation (4) may be converged after 3ns, but DA from Equation (5) still experiences significant shifts after this time (Figure 7c). Thel arger errors associated with the b3 system (ca. 0.2 kcal mol À1 ; Ta ble 4) appear to stem from the greater variation in V ses a valuess ampled,w ith ar ange in V ses a of 16 kcal mol À1 ,a sc ompared to values of 6a nd 11 kcal mol À1 for b1 and b2,r espectively.T his arises as af unctiono ft he steeper repulsive nature of the ses potentiala ppliedt ob3,s uch that small changes in the dihedral can lead to larger changes in energy.
The origin of the 0.70 AE 0.004 kcal mol À1 differencei nf ree energy for b3 butane in states 0a nd 1w hen using IT-TI is im- www.chemphyschem.org mediately apparent from the complete absence of dihedral transitions found in these simulations (Table 3a sw ell as Figures 3b,4ca nd 5c). Consequently, DA computedf rom IT-TI corresponds to the differencei ns tability of a t and g À conformer of butane. However,u nder the influence of the cooperative swarm of replica trajectories, the sesTI simulations sample the three rotamers of butane through each replica( Figures3ba nd 4c)a nd for each l (Figures 3a and 5c). Under the ses potential appliedt ob3,t he highest frequency of dihedralt ransitions for the butanem odelsi sf ound, with av alue of 35.2 ns À1 for f 0 (Table 3), permitting the sampling required to obtain estimates of free energy close to zero for the butane-to-butane transformation.

Discussion
Free energy calculations constitute valuablet ools in modelling chemicalp rocesses, for example, in predicting protonation state, solute partitioning betweeni mmiscible liquids and moleculara ssociation in the condensed phase. Free energy calculations, however,a re prone to an umber of potentials ources of error,p rincipally from the choice of model for the molecular Figure 3. Population of f 0 f 1 rotamerso fbutane during IT-TI or sesTIt ransformation in b1, b2 and b3 systems a) as af unction of l (combined replicas) and b) as afunctiono fr eplicasr 01 to r12 (for l = 0.01)a nd their sum ("all"). Abscissaisf 0 and ordinate is f 1 .Both axes range from À1808 to 1808.
ChemPhysChem 2015, 16,3233 -3241 www.chemphyschem.org system (e.g. force field, solvent modela nd treatment of electrostatic interactions) and from finite sampling. By using the butane-to-butane transformation,w ef ocus on the consequences of limited exploration of phase space. Here, errors arise from the omission of conformational regionst hat are important contributors to the differencei nf ree energy between two states. In the three butane modelsc onsidered, the t, g + and g À energy wells are low lying and contribute thermodynamically to an overall DA for ab utane-to-butane transition of zero. For b1,t he moderate energy barrier between t and g + /g À rotamers is adequately traversed by unbiased MD simulations of 3nsf or each l window.F or the 12 replicaI T-TI method, we obtain a DA estimate of À0.01 kcal mol À1 (Table 2). However, for the increasinglyk inetically distinct b2 and b3 models, the quality of sampling erodes and consequently the free energy estimates deviate to 0.21 and 0.70 kcal mol À1 ,r espectively.I n the latter case, entire rotamers are omitted. Accessing these important,b ut sometimesh idden, substates can be am ajor issue for free energy methodologies. [4,36] As an approacht or educing this source of error,w ee valuated aT I approachb ased on sesMD for conformational sampling. Multiple MD simulations of butanea re coupledt hrough their torsion angles by using as es potential with attractive and repulsive components [Eq. (3)].T he resultant sesTI calculations of DA for b1, b2 and b3 provideg ood sampling of all three butaner otamers and estimates close to zero ( Table 2). More frequentb arrierc rossing is also observed ( Table 3) as swarm replicast ransition between wells and stimulate crossings in neighbouring replicas. For b1,s esTI convergencea ppearsi mprovedo ver TI or IT-TI calculations, such that shorter l window simulation times can be employed. For b2 and b3,I T-TI calculations converge more quickly than sesTI, but to an incorrect pseudo-converged value, thus providing precision,t hat is, al ower statistical error associated with the trapped states, but not accuracy,o wing to significant systematic error.C onsequently,c orrection of the IT-TI estimates for b2 would require the application of much longer unbiased simulations for each l window,w hich appears beyond reach for the energy landscape of b3.
The improved sampling of butane rotamers in states 0a nd 1b yu sing sesTI is evident across l and the replicas (Figure 3). For all 13 l windows,e ach rotamerw ithin the space of f 0 and f 1 are sampled for b1, b2 and b3 by using sesTI, contrasting with sporadic transitions for b2 and no transitions for b3 (Figure 3a). This good coverage in sampling over l is the result of slightly differing f 0 f 1 distributions (averaged over l)f or each of the sesMD replicasr 01 to r12 (Figure 3b). The resulting aggregate distributions ("all" in Figure 3) show ab roader sampling arounde ach rotamer compared to IT-TI, sampling the higher-energy sideso ft he wells. This would be particularly important for situations where the minimaa re located atd ifferent geometries in states 0and 1.
For sesTI, we adopted two differentr e-weightings chemes, according to Equation (4) or Equation (5). In the first scheme, www.chemphyschem.org derived from an assumption of weak coupling between replicas, the replicas are re-weighted accordingt ot he biasingi nfluence of the ses potential, V ses a .Equation (4) re-weights on areplica basis, such that the dominance of any single replica to the computed average is effectively minimised. In the second scheme,d erived from the work of Malevanets and Wodak, [37] the further assumption is that of ergodicity in the replica trajectories,s uch that the time average over ar eplica will converge to the average over replicas. Resultantly, the average computed by using Equation (5) can accentuate the dominance of individual replica contributionsto @V MM ðlÞ @l l .We find that the two weighting schemes are in general agreement with each other,e xcept for the case of the highest energy barrier system, b3.Here, it is possible that the strong potentialpolarises the travel of certain replicas (Figure 3b), leadingt o as malld ivergence in the overall DA estimate of 0.11kcal mol À1 (Table 2). Indeed, this also underlies the longer convergence required for b3 through this re-weightings cheme( Figure7c).
Other free energy calculations with enhanceds ampling MD methods have been applied in addressing the issue of pseudoergodicity and constitute promising alternatives;t hese include well-tempered metadynamics, [38] replicae xchange with solute tempering (REST2) [39] and windowed aMD in aHamiltonian replica exchange framework (w-REXAMD). [40] Several biased MD approaches such as aMD use exponential re-weighting;f or systems with largeb iasinge nergies, broad distributions of these energies lead to high energy terms with sizeable exponential weights that dominate the free energy estimates. Conversely, the weights of low energy terms are often lost in the limitations of numerical precision.T here is indeed evidenceo fthis energetic noise in the application of the stronger ses potential for b3 here, which resulted in the broadest range in V ses a of the three systems. For applicationsw here more dihedrals are enhanced by the ses potential, the largerb iasinge nergy could potentially furtheri ncreaset he spread of V ses a values and the accompanying uncertaintyi nf ree energy estimates. Clearly, furtherw ork is required to assess the optimal ses parameters and limiting system size for the recovery of accurate free energy profiles. In this regard, we note that the use of approximations to the exponential term such as cumulant expansion [41] have shown utility in reducing noise in re-weighting, [42] such that an arrow distribution is maintained, even for the enhancement of ag reater number of degrees of freedom.
As ar eplica-based approach, sesTI affords as traightforward couplingo fp otentials, albeit with aj udicious choice of param- Table 3. Frequencyo ft ransitions of C-C-C-C dihedrals f 0 and f 1 in b1, b2 and b3 from IT-TI and sesTI calculations,a veraged over replicaa nd l. Standard deviations in parentheses.

Model
Method www.chemphyschem.org eters, and circumventst he energy overlap requirement of replica-exchange methods. In common with metadynamics, the sesTI approachr equires ac hoice of coordinates to sample, in its currentf orm as electiono fd ihedrala ngles. In principle, detailed knowledge of hidden, thermodynamically relevant conformations is not required ap riori, butinsteadc an be explored by using the swarm of coupled trajectories.

Conclusions
We have described af ree energy simulationa pproachb ased on as warmo fc oupled replicas to improvet he underlying sampling of kinetically distinct states. Computational free energy changes for this transition using dual topology TI and IT-TI increasingly deviated from zero with increasing barrier height of intramolecular rotation. Alternatively,d ualt opology sesTI calculations applied enhanced sampling to the intramolecular dihedral of the reference and target states of butane, and led to computed free energy differenceso fz ero for barrier heightsu pt o6kcal mol À1 .T he sesMD simulations underlying these improved free energy change estimates displayed increased frequency of transitions between wells and greater coverage of phase space,a ss warm replicas interacted to drive each other across energy barriers. The sesTI approach, therefore, shows potential in quantifying free energy differencesi n systemsw here the free energy change may depend on multiple distincts ubstates separated by non-trivial energy barriers. Figure 6. Convergence of free energy difference estimates as af unction of simulation length of 12 independent replica TI calculations for systems a) b1,b)b2 and c) b3. Figure 7. Convergence of free energy differencee stimates as af unction of simulation length for butane-to-butane transitionso fa)b1,b )b2 and c) b3 systems from IT-TI (blue)a nd sesTI by using Equation(4) (red) and Equation (5) (green).