the free - energy change upon binding is the fundamental thermodynamic quantity to evaluate inhibitor affinity for a target protein .
reliable free - energy changes can be estimated by computer simulations via thermodynamic - integration ( ti ) methods . in practice ,
such calculations are highly accurate for small compounds within the force field and model resolution employed . in principle
however , ti approaches require a sufficient sampling of the phase - space regions where the hamiltonians corresponding to two states of the system differ significantly .
therefore , the practical use of ti - based approaches in the context of macromolecular processes is still rather limited .
finite sampling problems for a given equilibrium thermodynamic state can be alleviated by multiple independent simulations .
this enhances phase - space sampling and allows distribution of the computation into a number of independent runs , which is particularly appealing in consideration of the rapid and steady increase of computational power in the form of multiple cpu clusters vs single cpu supercomputers ( e.g. , http://www.sdsc.edu ; http://www.nccs.gov ; http://www.bsc.es ) . here , we present the independent - trajectories thermodynamic - integration ( it - ti ) approach to calculate free - energy changes for ( bio)molecular systems .
it - ti employs multiple , independent ti calculations to calculate a free - energy change of interest , while incorporating both soft - core potentials and ligand translational restraints to effectively improve the extent of phase - space accessed .
our results show that it - ti allows significantly increased accuracy compared with standard ti . using it - ti in the context of proteinligand binding and macromolecular
this is the case for the h5n1 avian influenza neuraminidase receptor studied in this work ( figures 1 and 2 ) .
( a ) n1-pvr model structure and scaffold structures containing important pvr - binding residues ( green , -sheet 406 ; red , -turn 277 ; orange , -turn 222 ; pink , loop 119 ; blue , loop 150 ) .
time series of the c atom - positional root - mean - square deviation ( rmsd ) of each n1-pvr binding site loop or secondary structure element from the x - ray structure .
the avian influenza virus type a , particularly its h5n1 form , is becoming a worldwide pandemic threat due to its high virulence and lethality in birds , rapidly expanding host reservoir , and exceptionally elevated mutation rate ( http://www.who.int/csr/disease/avian_influenza ) .
extraordinary research efforts are devoted to understanding the molecular basis of inhibitor susceptibility to avian influenza viral enzyme neuraminidase ( na ) mutations , particularly for the lethal and drug - resistant group 1 na enzymes that include h5n1 .
the inhibitor peramivir ( pvr , also known as bcx-1812 or rwj-270201 ; developed by biocryst pharmaceuticals , birmingham , al ; see scheme 1 ) is demonstrated to be active in vitro and in vivo against both group 1 and 2 viral na .
therefore , pvr constitutes a promising candidate for further drug - design research.(24 ) the it - ti free - energy changes due to pvr hydration ghydr(l ) and n1-pvr binding gbind(l ) are shown as well as pvr alchemical modification ghydr(l * ) and gbind(l * ) .
all values are given in kj mol with corresponding g uncertainties between parentheses . in this paper
, we explore the changes of conformational dynamics and hydration of pvr upon binding to avian influenza virus h5n1 na ( figure 1 ) .
we perform it - ti calculations that yield an accurate estimate for the n1-pvr free energy of binding , within 1 kj mol of experiment . then , we investigate n1-pvr binding determinants and quantify their thermodynamic role in the binding process through it - ti alchemical modifications of selected pvr moieties .
this work represents a first step in the computer - based development of a putative novel class of n1 inhibitors from accurate free - energy calculations .
we anticipate that it - ti will allow , in general , the estimation of relative free energies for systems of increasing size , with improved reliability , by employing large - scale distributed computing.(25 )
the initial coordinates for the n1 neuraminidase monomer bound to the pvr inhibitor ( n1-pvr ) were taken from the x - ray crystal structure(26 ) of n1 bound to oseltamivir ( pdb i d : 2hu4 ; chain a ) , because no n1-pvr structure has been deposited to date .
atom positional coordinates for pvr were taken from the corresponding n8-pvr structure(26 ) ( pdb i d : 2htu ; chain a ) and superimposed onto 2hu4 using the protein backbone c atoms ; this results in superimposition of the oseltamivir and pvr ring atoms .
the n1-pvr complex was solvated in ( pre - equilibrated ) cubic boxes large enough ( 8.3 nm ) to avoid any interactions between mirror images under rectangular periodic boundary conditions .
three randomly chosen water molecules ( minimum ionion distances of 1.0 nm ) were replaced with na ions to neutralize the system .
initial configurations for the water ( wt - pvr ) and vacuum ( vc - pvr ) reference states were defined using the same pvr coordinates . for a summary of n1-pvr , wt - pvr , and vc - pvr simulated systems ,
all simulations were performed using the gromos05 software for biomolecular simulation(27 ) and the gromos force field(28 ) ( 45a3 parameter set(7 ) ) .
gromos pvr force - field parameters were derived from existing building blocks ( supporting information , table s1 ) .
the gromos compatible spc water model(31 ) and previously reported spc - water compatible parameters for ions(32 ) were employed . for ligand simulations in vacuo ( vc - pvr )
the corresponding 45b3 parameter set was employed . for n1-pvr , a first steepest - descent energy minimization ( em ) was performed to relax solvent and ions , while protein atom positions were restrained by using a harmonic potential ( force constant k = 2.5  10 kj mol nm ) . a second em run without restraints eliminated any residual strain .
all em runs were extended until the energy change per step became < 0.5 kj mol .
the system was then brought to the reference temperature ( t = 300 k ) in six consecutive md periods of 25 ps ( 50 k increments ) . during the heating of n1-pvr ,
protein atom positions were restrained with a harmonic potential , using a k from 10 to 0 ( decreased in steps of 2.5  10 kj mol nm ) .
in addition , four independent md runs were initialized by reassigning random velocities from maxwellboltzmann distributions at 5 k. all five independent trajectories were extended ( at least 2 ns ) to reach equilibration of the separate system hamiltonian components .
independent trajectories for the wt - pvr and vc - pvr systems were similarly prepared .
one md run for each system was extended for 20 ns and used for conformational analysis ( see below ) .
newton s equations of motion were integrated using the leapfrog algorithm(33 ) with a 2-fs time step .
the shake algorithm(34 ) was applied to constrain all bond lengths ( relative geometric tolerance of 10 ) .
all simulations were carried out in the n , p , t ensemble ( reference pressure 1 atm ) by separately coupling the temperature of solute and solvent degrees of freedom to a 300 k heat bath(35 ) ( relaxation time 0.1 ps ) and by coupling the pressure ( estimated based on an atomic virial ) to a pressure bath(35 ) via isotropic coordinate scaling [ relaxation time 0.5 ps ; isothermal compressibility 4.574  10 ( kj mol nm ) ] .
nonbonded interactions in the range 0.00.8 nm were recalculated every time step and in the range 0.81.4 nm every five time steps and truncated at 1.4 nm .
a reaction - field correction was applied to account for the neglected interactions beyond 1.4 nm,(36 ) using a relative dielectric permittivity of 61 for the spc water model.(37 ) a fast grid - based pairlist - construction algorithm(38 ) was employed ( cell - mask edge of 0.4 nm ; atomic - level cutoff ) as implemented in the gromos05 md++ module.(27 ) trajectory snapshots were extracted every 2 ps from the 20-ns simulations .
structural fitting was performed by ( i ) superimposing solute centers of mass ( to remove overall translation ) and ( ii ) performing an atom - positional least - squares fitting procedure(39 ) ( to remove overall rotation ) using n1 c atoms or all pvr atoms .
transient n1-pvr interactions identified as important binding motifs were monitored using the gromos++ analysis software.(27 ) hydrogen bonds were defined to have a maximum hydrogen - acceptor distance of 0.3 nm and a minimum donor - hydrogen - acceptor angle of 125. an extended hydrogen bond criterion was used ( 0.35 nm ; 120 ) to capture additional relevant interactions .
salt bridges and hydrophobic contacts were considered formed for atom pair distances < 0.45 nm . secondary structure elements were defined by the following n1 residue sequences : loop 119 , v116p120 ; loop 150 , t148s153 ; loop 277 , i222e227 ; -turn 277 , h274c278 ; -turn 292 , v290n294 ; loop 347 , g345k350 ; loop 371 , s368g373 ; -sheet 406 , s404g408 ; loop 430 , r430w438 . the free - energy change between two states a and b can be estimated by thermodynamic integration ( ti ) as(40 ) 
 where h( ) denotes the system hamiltonian from a single trajectory as a function of the coupling parameter  and  ...  denotes ensemble averaging at a given  value . in it - ti ,
hi( ) is the system hamiltonian for the ith independent trajectory , and the mean free - energy change gab reads 

 where the integration runs over n independent trajectories . in principle , under the assumptions of ( i ) infinitely
long trajectories and ( ii ) a fully accessible system phase space , eq 1 will provide an estimate of the free - energy change between two states a and b which is identical to that provided by eq 2 .
this follows in the limits of validity of the ergodic hypothesis . in practice , however , due to the fact that ( i ) only finite simulation times can be achieved and ( ii ) the phase space of a solvated macromolecule is far from being fully accessible ( i.e. , its corresponding free - energy landscape is a very rough and frustrated surface at standard / physiological conditions ) , eqs 1 and 2 provide significantly different free energy estimates ( see ) .
it - ti overcomes this practical limitation by enhancing phase - space sampling of the thermodynamic systems of interest , therefore adding to the reliability and predictive power of free - energy calculations .
two types of thermodynamic perturbations ab were performed in this study , alternatively employing eq 1 or 2 : ( i ) from ligand l full potential (  = 0 ) to zero nonbonding interactions (  = 1 ) ; ( ii ) from the l full potential (  = 0 ) to that of a chemically modified ligand l * (  = 1 ) ; see scheme 2 . in both cases , soft - core interaction potentials(16 ) were used for l atoms involved in the perturbation ( slj = 0.5 and sc = 0.5 ) to avoid singularities and to enhance phase - space sampling .
the binding free energy gbind(l ) for ligand l can be estimated by an annihilation perturbation (  = 0   = 1 ) of the nonbonded ligand interactions in the protein and water - solvated thermodynamic states .
the impact on the binding free energy , gbind(l * ) , for a ligand modification l * can be estimated by corresponding perturbations ( l  l * ) in both thermodynamic states .
two alternative statistical procedures were employed to evaluate the uncertainty  for gab or gab free - energy estimates .
first , a simulation standard error sim(t ) of the time - varying hamiltonian derivative at a given  can be calculated as 

 with t being the total number of block averages(41 ) throughout the single ith trajectory or all n concatenated independent trajectories .
( ht()/) denotes the hamiltonian derivative , block - averaged at time t , and ht()/ is the ensemble average over the entire simulation time at a given . as an example , sim(t ) uncertainties are reported as error bars for ht()/ vs  in figure 3 ( solid black curve ) .
then , a corresponding free - energy uncertainty can be obtained as 
  h()/ values and corresponding uncertainties sim(t ) vs . ( a ) annihilation of pvr in n1 binding site ( gn1(l ) ; scheme 2a ) .
( b ) annihilation of pvr in water ( gwt(l ) ; scheme 2a ) .
inset panels highlight  regions where it - ti averages outperform standard individual ti calculations .
this follows from the standard assumption that ( ht()/) values are statistically uncorrelated along the time over different values of the coupling parameter . however , the gi uncertainty includes the physically based fluctuations of ( ht()/) , though corresponding noise is typically reduced by block - averaging.(41 ) therefore , despite its wide use in the literature , gi is a questionable measure of uncertainty for a free - energy change of interest .
for example , considering that overlap of phase space at neighboring  values is a requirement for smooth ht()/ vs  curves ( figure 3 ) , one could claim that ( ht()/) time series are statistically correlated .
nonetheless , the abovementioned uncertainty defined in eq 4 is representing the lowest possible uncertainty for a free - energy - change estimate from standard ti .
thus , it seems the fairest choice for this study comparing ti vs it - ti results .
second , for it - ti , a statistics - based uncertainty g on a given free - energy change gab from eq 2 can be calculated as the standard deviation from the mean ( standard error ) of the n gi results 

 where g is the standard deviation of the free - energy change over the n it - ti trajectories employed .
importantly , g has a clear statistical validity,(42 ) because of its explicit dependence on the repeated independent estimates .
similarly , for a general overall free - energy change gab , calculated as the difference between two free - energy changes gb and ga , a corresponding uncertainty can be obtained by propagating the respective uncertainties as(42 ) then , the relative uncertainty for a given free - energy change a  b reads 
  in this study , it - ti runs were extended to obtain sufficiently smooth curves of hi()/ vs  ( figure 3 ) .
it - ti trajectories were independently equilibrated ( 0.5 ns ) for each of the 26  points ( from five initial equilibrated  = 0 configurations ) , followed by independent sampling periods ( 0.5 ns ) used for free - energy estimation . increased
sampling ( up to 2.5 ns ) times were required in the ranges 0.12    0.24 and 0.76    0.92 .
only n1tail1 and n1tail2 modification perturbations had larger g(% ) values ( up to 52% ) due to the corresponding small gab values ( supporting information , table s3 ) . for gn1(l ) ( scheme 2a ) , the potential u(rl ) = /2k(rl  r0 ) was applied to harmonically restrain ligand translation and ensure its sampling of a finite phase - space volume v. an optimal k value of 246.5 kj mol nm was estimated from the ensemble - averaged l root - mean - square deviation ( rmsd ) by applying the equipartition theorem as(18 ) 
 where rl is the position of l during 20-ns of unrestrained simulation at  = 0 , r0 is the initial position of l in the pre - equilibrated starting configuration , p(rl ) is the 3-d positional probability distribution of l , r is the molar gas constant , and t = 300 k. in this study , rl and r0 are the instantaneous and initial positions of the pvr c4 ring carbon .
the phase - space volume v can be defined as(18 ) 
 where  = 1/rt .
thus , the correction(17 ) to the restraining potential bias for gn1(l ) in scheme 2 reads 

figure 1 shows the six - bladed -propeller structure of n1 neuraminidase bound to pvr and important active site residues , monitored throughout 20 ns of n1-pvr md simulation . to investigate n1-pvr receptor flexibility and conformational changes , the c atom - positional rmsd from the x - ray crystal structure was monitored for the secondary structure elements forming the n1-pvr binding site ( figure 2 ) .
loop 150 deviates the most ( up to 0.7 nm ) to sample stable open - loop configurations ( 0.5 units ) throughout the md simulation ( figures 1 and 2 , blue ) .
this observation is similar to that recently reported for both apo and oseltamivir - bound n1 neuraminidase ( up to 0.6 nm).(43 ) additionally , we find that loop 119 demonstrates significant flexibility and samples conformations with rmsd of 0.4 and 0.3 nm in the n1-pvr complex ( figures 1 and 2 , pink ) .
the dynamics of loops 150 and 119 in water indicate a significant relaxation from the crystal - packed conformation captured in x - ray experiments .
the -sheet 406 and -turns 277 and 222 show comparatively lower rmsd deviations ( 0.10.3 nm ) and smaller fluctuations on the 20-ns time scale .
a dominant multicenter hydrogen bond between n1 e277 and the pvr guanidinium group is stable for the entire simulation time , with pvr nr3 and nr1 atoms alternating as hydrogen donors to e277 carboxyl oxygen atoms oe1 and oe2 ( 5069% occurrences ) .
e277 also transiently interacts with the pvr methyl acetamide polar hydrogen ( hp ; 57% ) . yet
another hydrogen bond ( 54% ) is observed between the pvr carboxyl oxygen ( od1 ) and the n1 y406 hydroxyl ( oh ) .
interestingly , the y406 residue is homologous to other key catalytic tyrosine residues found in the avian influenza virus family of glycosidases ( carbohydrate active enzymes database ; http://www.cazy.org/ ) .
thus , targeting y406 could be important for drug design , as suggested by a computational solvent mapping analysis.(44 ) an extended hydrogen bond between n1 r152 on loop 150 and pvr ( oh group ) occurs for 44% of the simulation time ( table 2 and figure 1 ) and contributes to the stability of open loop 150 ensemble of configurations
oe refers to glutamate carboxyl , oh to tyrosine hydroxyl , and nh2h22 to arginine guanidinium .
all n1-pvr intermolecular hydrogen bonds occurring > 5% are shown . using extended hydrogen - bonding criterion .
n1 protein atoms that maintain strong interactions with either ( i ) charged ligand moieties ( salt bridges ) or ( ii ) hydrophobic groups ( hydrophobic packing ) display distance probability distributions with one sharp , high - intensity peak .
those residues experiencing varied interactions have broader distributions and/or multiple peaks . throughout the 20-ns simulation , the pvr guanidinium group ( cr atom , scheme 1 )
maintains a well - defined salt - bridging interaction with the n1 e277 carboxyl group ( 0.38 nm average distance ) .
the broad distribution for the pvr crn1 e119 distance ( 0.52 nm ) and the doubly peaked distribution for the pvr crn1 e227 distance ( 0.47 and 0.65 nm for first and second peak ) correspond to transient salt bridges formed by these n1 glutamate residues to the pvr guanidinium group .
conformational sampling of pvr moieties in the n1 neuraminidase binding site from 20-ns n1-pvr molecular dynamics simulation .
distance pairs are labeled according to the atoms monitored : pvr cr to n1 e277 ( blue ) , e119 ( green ) , and e227 ( red ) carboxyl group carbon ; pvr cd to n1 r118 ( pink ) , r292 ( yellow ) , and r371 ( violet ) ; pvr c11 to n1 i222 side chain carbon ( black ) and w178 aromatic ring carbon ( cyan ) .
stable long - range interactions between pvr and the arginine triad of r118 , r371 , and r292 residues ( figure 1 ) are observed along our 20-ns n1-pvr simulation , with distribution peaks at 0.68 , 0.8 , and 0.96 nm , respectively ( figure 4 ) .
this triad of positively charged residues also contributes to binding via the conserved ligand carboxyl group in both the natural ligand , sialic acid , and other synthetic neuraminidases inhibitors ( e.g. , dana , oseltamivir , and zanamivir ) . our results support the fundamental role of r118 and r371 in n1-pvr binding .
the distal position of r292 indicates less reliance of pvr on this residue and could account for the retained pvr affinity for the n2 resistance mutation r292k .
the pvr aliphatic tail ( scheme 1 ) has been designed(45 ) to fill a small hydrophobic subpocket , comprised of w178 and i222 residues , which is conserved among both group 1 and 2 na ( figure 1 ) .
our md simulation confirms that pvr tail atoms interact with the branched i222 side chain ( c11 peak at 0.40 nm average distance ; figure 4 ) .
the same pvr tail atoms are also stably close to the w178 aromatic carbons ( c11 peak at 0.42 nm ; figure 4 ) .
these results highlight the occurrence of important hydrophobic - packing interactions in the n1 subpocket .
we note that the majority of the conserved residues described above have been suggested to also participate in key interactions between sialic acid and inhibitors dana , zanamivir , and oseltamivir , as well as pvr , in both group 1 and 2 na receptors .
figure 5 shows the radial distribution functions ( rdf ) for water oxygen atoms from wt - pvr and n1-pvr simulations .
the pvr guanidinium group ( cr atom ; scheme 1 ) undergoes 0.5 units decrease in its first peak intensity upon n1-pvr binding .
we can explain this desolvation effect by considering the tight interaction of this bulky , positively charged pvr group with the negatively charged n1 e277 carboxyl ( figure 4 and table 2 ) .
desolvation upon n1 binding is also observed for the pvr aliphatic tail by the intensity decrease of its first solvation shell peaks ( c10 and c11 atoms , 0.50.9 units ; figure 5 ) .
this is consistent with the formation of more favorable interactions between pvr and n1 residues i222 and w178 .
ensemble - averaged solvation of pvr moieties when free in solution ( left panel ) or bound to the n1 protein binding site ( right panel ) .
radial distribution functions of the water oxygen atoms from 20-ns molecular dynamics simulations are shown centered on pvr atoms cr ( blue ) , c10 ( black solid ) , c11 ( black dashed ) , cd ( red ) , o5 ( green ) , and np ( cyan ) .
a different hydration behavior can be noticed for the pvr carboxyl group ( cd atom ) , with a similar strong intensity for its first peak at 0.35 nm in both wt - pvr and n1-pvr simulations .
thus , the pvr carboxyl group is still solvated in the n1 binding site by dynamic water molecules on an ensemble averaged basis .
the presence of water molecules in the n1-pvr binding site and the lack of persistent ligandsolvent hydrogen bonds confirm this point .
a pvr carboxyl groupn1 y406 hydroxyl hydrogen bond is transiently formed ( 54% occurrence ; table 2 ) , allowing this pvr moiety to still repeatedly interact in the n1-pvr binding site with water molecules .
water exchange in charged protein cavities(49 ) and water - mediated interactions in flexible carbohydrateprotein binding(50 ) have been previously reported .
the pvr hydroxyl group displays a first peak with reduced intensity when n1-bound ( o5 atom , 0.5 units ; figure 5 ) .
this moiety forms a competing hydrogen bond to n1 r152 ( 44% ; table 2 ) , yet water molecules are maintained in the first solvation shell upon binding .
the pvr acetamide group ( np atom ) has limited solvent accessibility due to the adjacent hydrophobic tail .
its solvation is further decreased ( 0.5 units ) in the bound state upon formation of a hydrogen bond to n1 e277 ( 57% occurrence ; table 2 ) .
overall , the interactions of the pvr guanidinium group with n1 e277 , e227 , and e119 ; the pvr carboxyl group with catalytic n1 y406 ; and the pvr aliphatic tail with n1 w178 and i222 in the hydrophobic subpocket appear most relevant to drive n1-pvr binding based on conformational and hydration analyses .
the free - energy change upon n1-pvr binding was estimated using the it - ti method and compared with standard ti values from single trajectories , as well as with experiment .
examples of hi()/ vs  curves are shown for standard ti and the improved it - ti calculations of the n1 protein and water reference states ( figure 3 ) .
the it - ti curves are smoother than those obtained from individual standard ti runs , because of increased sampling and improved overall statistics obtained through n = 5 independent ensembles .
their integration ( eq 2 ) provides gwt(l ) and gn1(l ) values to estimate gbind ( supporting information , table s3 ) .
table 3 reports the it - ti results and their comparison with the available experimental data .
our it - ti gbind estimate ( 61.1  5.4 kj mol ) matches the gbindexp value derived from multiple ic50 measurements ( 62.2  1.8 kj mol ; see supporting information , table s4 ) . a free - energy difference of 1.1
kj mol ( i.e. , 2% relative difference ) has no statistical significance within the above uncertainties .
remarkably , such an it - ti prediction of the experimental value relies on individual gbind(n ) estimates that span a rather large free - energy range ( see table 3 ) .
in fact , these standard ti estimates are at variance with the independent calculation performed , ranging from a substantial underestimation ( i = 4 ; 19% relative difference ) to a substantial overestimation ( i = 5 ; 29% relative difference ) of the gbindexp value .
only one of the individual standard ti results is in good agreement with the gbindexp value ( 63.7  11.0 kj mol , i = 1 ; 2% relative difference ) .
a significantly different n1-pvr binding free - energy estimate of 1180.9  31.8 kj mol has been reported based on mm - pbsa calculations.(48 ) the average gbindexp value from repeated experiments can be compared with the it - ti gbind estimate .
derived using data in refs ( 21 ) , ( 51 ) , and ( 52 ) ( see supporting information , table s4 ) . equation 5 , with propagated uncertainties as in eq 6 .
the it - ti free - energy estimate also has a lower g(% ) relative uncertainty compared to the gi(% ) from standard ti , i.e. 9% of the calculated free - energy difference ( table 3 ) .
gi(% ) values associated with independent gbind(n ) estimates are larger , ranging between 12% ( i = 4 ) and 19% ( i = 5 ) .
we stress that the uncertainties gi typically evaluated for single standard ti trajectories are not statistical indicators of the gbind accuracy ( see , eq 4 ) .
instead , repeated it - ti runs allow calculation of the more representative free - energy uncertainty g ( see , eq .
5 ) . in addition , for this study , the latter is directly comparable to the standard error of gbindexp determined from n experimental values ( supporting information , table s4 ) .
overall , our results underscore the improved predictive power of it - ti vs standard ti , due to the increased statistical reliability
. the large deviations observed among standard ti estimates can be explained , in part , by the flexibility of the n1 binding site .
loops 119 and 150 demonstrate heterogeneous conformational sampling among different  regions ; their rmsd from the initial equilibrated structure reach values up to 0.4 and 0.7 nm , respectively ( data not shown ) .
this is consistent with both the dynamic loop behavior from the longer 20-ns md simulation ( figure 2 ) and with the large gi values for standard ti estimates ( table 3 ) .
we conclude that the it - ti method significantly aids sampling of thermodynamic macrostates for flexible receptors by ensemble averaging of independent trajectories . n1-pvr
binding determinants were also investigated using it - ti free - energy changes upon computer alchemical modifications and their underlying thermodynamic cycle ( scheme 2b ) .
scheme 1 summarizes the corresponding gbind values together with the pvr free energy of binding gbind .
a positive or negative value of gbind indicates thermodynamically unfavorable or favorable alchemical modifications of the ligand l. neutralizing coo and nr3 charges has large but opposite effects on ligand binding ( gbind of + 55.1  3.1 and 79.7  4.2 kj mol , respectively ) .
the tail1 and tail2 modifications both have small , favorable impacts on ligand binding ( gbind of 5.8  2.4 and 1.5  2.0 kj mol , respectively ) . to understand these results , one must look at the effects on receptorligand interactions and ligand hydration free - energy changes ( ghydr ) upon alchemical perturbation . in other words ,
for example , a positive gbind value may be driven by ( i ) unfavorable ( enthalpic or entropic ) n1ligand interactions , ( ii ) a more favorable ghydr , or ( iii ) a thermodynamically unfavorable combination of the previous effects .
similarly , a negative gbind value may be driven by ( i ) favorable ( enthalpic or entropic ) n1ligand interactions , ( ii ) a more unfavorable ghydr , or ( iii ) a thermodynamically favorable combination of the previous effects . in this section ,
we address the impact of different it - ti modification perturbations on n1ligand binding ; in the next section , we consider the ligand hydration free energy . both are needed to fully describe a given gbind binding free - energy change . throughout the  = 0   = 1 n1coo modification perturbation , r292 and r371 residues move on average apart from the ligand scaffold .
the closest arginine , r118 , reduces its average distance to the ligand cd atom ( supporting information , figure s1a ) .
the important ligand carboxyly406 interaction ( table 2 ) is partially disrupted , while guanidinium interactions with e119 , e227 , and e277 are maintained . during the  = 0   = 1 n1nr3 modification ,
n1 r371 samples a more stable conformation close to the ligand cd atom , as revealed by a sharper distance distribution peak ( cf .
, the ligand guanidinium loses its favorable electrostatic interaction with n1 e277 , which shifts away from the perturbed moiety .
however , the ligand acetamide group is pushed closer to residue r156 , which forms a hydrogen bond with atom op ( data not shown ) .
this residue is not observed to closely interact with the unperturbed pvr molecule . following the  = 0   = 1 tail1 perturbation , i222 and w178 distance distributions for the modified pvr aliphatic tail ( c8 and c9 atoms ) transition to sharper and fewer peaks ( cf .
 = 0 vs  = 1 ; supporting information , figure s1c , d ) . on the other hand , the  = 0   = 1 tail2 perturbation has a limited effect , as the distributions of tail atoms c8 and c10 with i222 and w178 residues remain predominantly broad ( 0.40.9 nm ; supporting information , figure s1e , f ) .
scheme 1 summarizes the hydration free energy for pvr , ghydr , and the changes of this quantity for its alchemical modifications , ghydr .
the corresponding binding free - energy , gbind , and changes of this quantity , gbind , are also reported .
binding free energies ( gbind , gray solid bars ) and corresponding hydration free energies ( ghydr , diagonal lined bars ) for all modification perturbations in this study .
vertical bars display the corresponding uncertainties g. corresponding free - energy values and uncertainties can be found in scheme 1 and in supporting information , table s3 .
we estimate an it - ti ghydr value of 107.8  1.2 kj mol for the pvr molecule .
we note that this value is in qualitative agreement with hydration free energies of a large variety of compounds used to calibrate the force field used in this study .
the value of 820.5  40.2 kj mol reported based on mm - pbsa calculations significantly overestimates the favorable thermodynamic effect of pvr hydration.(48 ) the coo perturbation gives a positive ghydr value of 15.7  5.1 kj mol .
when deleting its charge , water structure is reduced (  = 0 vs  = 1 ; wt - coo ; supporting information , figure s2 ) .
ligand hydration around this moiety is maintained in the n1 binding site , yet its first solvation shell is displaced ( cf .
 = 0 vs  = 1 ; n1-coo ; supporting information , figure s2 ) .
a sufficiently unfavorable ligand hydration free energy in the unbound state would drive hydrophobic binding .
however , this effect is not large enough to compensate for the loss of favorable n1ligand interactions ( see the ) .
the coo charge perturbation has the largest unfavorable impact in the n1-bound state , leading to a gbind value of 55.1  3.1 kj mol .
a different thermodynamic compensation occurs between the bound and unbound states for the nr3 charge perturbation , with a large unfavorable ghydr value of + 124.9  1.8 kj mol .
the single solvation peak of the pvr guanidinium group in water shifts to larger distances ( cf .
 = 0 vs  = 1 ; wt - nr3 ; supporting information , figure s2 ) due to hydrophobic desolvation of this bulky charge group .
the limited hydration of the charged pvr guanidinium in the n1 active site is almost unaffected by the perturbation ( cf .
 = 0 vs  = 1 ; n1-nr3 ; supporting information , figure s2 ) .
the nr3 charge perturbation has the largest unfavorable impact on the unbound state , overcompensating for the loss of n1ligand favorable interactions ( see ) .
this leads to a net gbind change of 79.7  4.2 kj mol , significantly more favorable than the pvr gbind value .
the experimentally observed , improved binding of inhibitors to n2 neuraminidase by hydrophobic substitution at the pvr guanidinium position supports these results.(53 ) in water , tail1 and tail2 modifications increase the solvation around the unperturbed atoms ( cf .
 = 0 vs  = 1 ; supporting information , figure s2 ) and correlate to the ghydr values of 1.4  0.6 and 7.3  0.3 kj mol , respectively . in the n1 binding site ,
a distinct decrease ( 0.5 units ) for tail1 solvation ( c9 atom ; supporting information , figure s2 ) agrees with the rearrangement of n1 hydrophobic residues i222 and w178 ( see above ) .
the opposite signs of ghydr offset the changes observed in the protein as well as entropic changes to the aliphatic tail , resulting in similar gbind values of 5.8  2.4 and 1.5  2.0 kj mol for tail1 and tail2 alchemical modifications ( figure 6 ) .
overall , these results emphasize the dominant electrostatic contribution to the free energy of n1 binding for pvr and its alchemically modified variants and suggest that future drug development may also be guided by conveniently tuning ligand flexibility and hydrophobicity .
the independent - trajectories thermodynamic - integration ( it - ti ) approach was presented .
it allows for estimation of improved free - energy changes for biomolecular systems based on multiple independent simulations .
our results underscore the improved predictive power of it - ti vs standard ti , due to the increased statistical reliability .
standard ti estimates from individual trajectories span a rather large free - energy estimate range , from a 19% underestimation to a 29% overestimation of the experimental reference value ( 62.2  1.8 kj mol ) .
remarkably , our it - ti binding free - energy estimate ( 61.1  5.4 kj mol ) is in excellent agreement , i.e. 2% relative difference .
a general formulation is proposed to evaluate corresponding it - ti free - energy uncertainties that rely on a statistical treatment of error analysis .
overall , it - ti seems particularly promising in the case of highly flexible protein receptors , ligands , and macromolecular binding partners in general . using 20-ns molecular dynamics simulation of the n1-pvr complex
the interactions of the pvr guanidinium group with n1 e277 , e227 , and e119 ; the pvr carboxyl group with catalytic n1 y406 ; and the pvr aliphatic tail with n1 w178 and i222 in the hydrophobic subpocket appear most relevant to drive n1-pvr binding , based on conformational and hydration analyses .
furthermore , it - ti was applied to explore the binding determinants of avian influenza n1 neuraminidase inhibition using alchemical modification of the pvr molecule .
charge annihilation of its carboxyl and guanidinium groups has the largest unfavorable impact in the n1-bound and unbound states , respectively .
alchemical modifications of the pvr aliphatic tail suggest that future drug development may also be guided by conveniently tuning ligand flexibility and hydrophobicity .
finally , this study allows us more general conclusions on free - energy calculations in the context of proteinligand binding .
the key to designing improved inhibitors for a given target relies on an accurate thermodynamic description of both ligand - bound and ligand - unbound receptor and ligand states .
consequently , we suggest that the most reliable and predictive free - energy calculations will likely rely on the use of explicit solvent simulations and md force fields based also on a direct and general parametrization of solvation thermodynamics .
we anticipate the application of the it - ti approach to develop improved and potent drugs to inhibit flexible macromolecular receptors .
 free - energy changes are essential physicochemical quantities for understanding most biochemical processes . yet , the application of accurate thermodynamic - integration ( ti ) computation to biological and macromolecular systems is limited by finite - sampling artifacts . in this paper 
 , we employ independent - trajectories thermodynamic - integration ( it - ti ) computation to estimate improved free - energy changes and their uncertainties for ( bio)molecular systems . 
 it - ti aids sampling statistics of the thermodynamic macrostates for flexible associating partners by ensemble averaging of multiple , independent simulation trajectories . 
 we study peramivir ( pvr ) inhibition of the h5n1 avian influenza virus neuraminidase flexible receptor ( n1 ) . 
 binding site loops 150 and 119 are highly mobile , as revealed by n1-pvr 20-ns molecular dynamics . due to such heterogeneous sampling , standard ti 
 binding free - energy estimates span a rather large free - energy range , from a 19% underestimation to a 29% overestimation of the experimental reference value ( 62.2  1.8 kj mol1 ) . 
 remarkably , our it - ti binding free - energy estimate ( 61.1  5.4 kj mol1 ) agrees with a 2% relative difference . in addition , it - ti runs provide a statistics - based free - energy uncertainty for the process of interest . using 800 ns of overall sampling , we investigate n1-pvr binding determinants by it - ti alchemical modifications of pvr moieties . 
 these results emphasize the dominant electrostatic contribution , particularly through the n1 e277pvr guanidinium interaction . 
 future drug development may be also guided by properly tuning ligand flexibility and hydrophobicity . 
 it - ti will allow estimation of relative free energies for systems of increasing size , with improved reliability by employing large - scale distributed computing . 
Introduction
Materials and Methods
Results and Discussion
Conclusion