the 
 binding of a drug - like molecule to a protein leads to displacement 
 of water molecules from the protein s binding pocket , and the 
 thermodynamics of this displacement process is thought to contribute 
 significantly to the overall thermodynamics of protein 
protein hydrogen bonds may incur a large energetic penalty , 
 whereas displacement of water from hydrophobic parts of the binding 
 pocket may help drive ligand - binding . intuitively , one may view different 
 parts of the protein s surface as imposing different surface 
 energies on the nearby water , with correspondingly different thermodynamic 
 consequences for water displacement by various ligands .
the 
 use of molecular distribution functions to analyze molecular dynamics ( md ) simulations has led to important 
 advances in the study of binding site water and its role in molecular 
 recognition ; parallel progress with the 3d rism approach also deserves mention but is not considered here .
key early contributions 
 include development of watermap ( schrdinger 
 llc ) , stow , and other approaches , which have provided new insight and shown promise as tools to help 
 discover small molecules that will bind a targeted binding pocket . 

such methods frequently define spherical sites , where water is present 
 at high density , to represent the distribution of water in the binding 
 site .
this hydration site approach ( hsa ) is motivated in part by the 
 practical consideration that , in regions where water is present at 
 lower density , it becomes more difficult to obtain converged values 
 of the local orientational entropy of water .
this is a simple consequence 
 of the lower number of water samples available from the simulation 
 in such low - density locations .
the hsa strategy of limiting attention 
 to hydration sites where water is present at high density maximizes 
 the chances for good numerical convergence of the orientational entropy . 

however , as previously discussed , the 
 regions in a binding site where water is present at high density can 
 have a complex shape , which is not easily represented by a collection 
 of spheres .
this limitation has been addressed in a grid - based 
 implementation 
 of inhomogeneous solvation theory ( ist ) , termed gist . instead of constructing hydration sites
, gist discretizes the smooth 
 distributions of water density and other properties onto a fine , three - dimensional 
 grid .
the problem of converging the local orientational entropy of 
 water is overcome through the use of a highly efficient nearest - neighbor 
 ( nn ) method , as opposed to histogram methods , which require more sampling 
 to reach adequate convergence .
gist can also take 
 advantage of the fact that regions of lower density contribute proportionately 
 less than regions of higher density regions to the overall orientational 
 entropy of the displaced water .
this density - weighting means that , 
 if one is interested in the integral of the orientational entropy 
 over a volume containing both high and low density regions , one can 
 converge the overall integral to an acceptable tolerance , so long 
 as the high - density regions are well converged .
alternatively , the 
 grid approach makes it straightforward to focus on regions where water 
 is present at high density , as done in hsa , without simplifying their 
 shapes . here
protein 
 system . in order to establish a clear basis for comparing methods , 
 we study coagulation factor xa ( fxa ) with a set of small molecule 
 inhibitors used in early studies of the watermap method ,
and we derive scoring functions based on both 
 gist and hsa methods . for this initial test of gist
s applicability 
 to protein  ligand modeling , we do not seek to establish a full - fledged 
 protein  ligand scoring function , suitable for virtual screening 
 or lead optimization .
instead , as previously done , we ask how well the gist treatment of hydration can capture 
 affinity differences between closely related congeneric pairs of ligands , 
 where differences in binding affinity that result from contributions 
 other than solvation such as configurational entropy and protein 
our results support the applicability of gist 
 and , in addition , provide an unexpected outcome regarding the role 
 of energetically versus entropically unfavorable water .
we ran explicit - water md simulations of fxa and 
 used both gist 
 and our local hsa implementation to extract information about the 
 structure and thermodynamics of the water in the binding site .
we 
 then considered the displacement of this binding site water by various 
 fxa inhibitors , whose binding site poses are known or could be inferred 
 from the known poses of very similar compounds .
candidate scoring 
 functions , based on the computed properties of the water displaced 
 by each ligand , were trained on subsets of the experimental affinity 
 data and then tested on separate sets , in order to assess the utility 
 of the hydration data to resolve the relative binding affinities of 
 pairs of congeneric ligands .
details of the computational methods 
 are presented in the following subsections . as previously 
 detailed
, gist uses a three - dimensional rectangular grid of cubic 
 voxels in the region of interest and processes the snapshots of an 
 md trajectory to compute the following thermodynamic quantities for 
 each voxel , k , centered at location rk :  ek , swnorm , the mean solute  water interaction energy of a 
 water molecule in voxel k.  ek , wwnorm , one - half the 
 mean interaction energy of the water in voxel k with 
 all other waters .
the factor of 1/2 , which was not included in the 
 definition of the corresponding water 
water interaction energy 
 in our initial presentation of gist , is 
 customary in liquid - state theory ; it allows the total energy of neat 
 water to be written as the sum of the individual water energies .
 tsk , swtrans , norm , the single - body ( one - water ) translational entropy of water in voxel k , relative to bulk , normalized to the mean number of waters 
 in the voxel .
 tsk , sworient , norm , the single - body ( one - water ) orientational 
 entropy of water in voxel k , relative to bulk , normalized 
 to the mean number of waters 
 in the voxel .
note that we previously used a somewhat different 
 notation for 
 these quantities ; for example , ek , swnorm was eswnorm(rk ) . here ,
a 20.5  20.5  22.5 
  grid was centered on the active site of fxa .
the grid spacing 
 of 0.5  provides voxels large enough to give statistically meaningful 
 data but small enough to still give a high resolution description 
 of the density distribution functions .
it is worth noting that the volume of each voxel is 33.5 times less 
 than that of a hydration site ( see below ) , as the latter represents 
 a sphere of radius 1  .
the gist lennard - jones and electrostatic 
 energies were computed from stored md frames , using the minimum image 
 convention and no cutoff , and the reference value of the bulk water 
the main gist 
 calculations presented here used 100  000 frames saved at 1 
 ps intervals during a 100 ns production md run , but shorter durations 
 were also examined , to study convergence , as detailed below .
when 
 a ligand binds a protein , it displaces water from the protein s 
 binding site .
if the displaced water was unfavorable relative to bulk , 
 then water displacement should make a favorable contribution to the 
 ligand s binding affinity . with this in mind
, we initially 
 looked for a correlation between measured ligand binding affinities 
 and regional hydration free energies ( eq 25 of ref ( 28 ) ) , where the region was 
 defined as those voxels covered by each ligand in a bound pose .
however , 
 finding little correlation ( data not shown ) , we conjectured that any 
 underlying correlation had been obscured by noise , due to sharp variations 
 in the hydration energies with even small changes in ligand position
. 
 this sharpness traces , at least in part , to our use of a single pose 
 for each ligand , and a restrained protein structure in the water simulations . 

we addressed this issue by constructing scoring functions which are 
 based on the gist data but are less sensitive to the details of local 
 water properties , due to the use of cutoffs in local water density , 
 energy , and entropy .
the use of cutoffs to construct a well - behaved 
 scoring function from local hydration data was first introduced in 
 the context of a hydration site model .
we tested three such scoring functions based on the gist hydration 
 data available from the grid described above . in all three cases
, 
 voxel k can contribute to a ligand s score 
 only if the voxel s center , rk , lies within the van der waals radius of any atom 
 of the ligand . for a given ligand , then , each voxel is assigned a 
 binary displacement indicator , dk , which equals 1 if the center of the voxel lies within the 
 van der waals radius of any ligand atom and 0 if it does not . the 
 van der waals radii
are drawn from the software package crystal maker 
 ( crystalmaker software ltd . ) , which in turn relies on bondi .
we also allowed the scoring function to focus 
 on voxels where water is present at high density by setting up an 
 additional binary indicator , gk , which is set to 1 if the water density in voxel k exceeds a cutoff gco , and 0 if it does 
 not .
this cutoff is one of the trained parameters , so it will be greater 
 than zero only if imposing a density cutoff actually improves the 
 accuracy of the scoring function . finally , we set up similar cutoffs 
 for the total energy and entropy , eknorm and tsk , swnorm , associated with each voxel k , and used these to define additional binary masks based 
 on energy and entropy thresholds .
thus , the binary mask , ek , equals 1 if eknorm exceeds the cutoff , eco , and 0 otherwise ; 
 and the binary mask , sk , equals 1 if tsk , swnorm exceeds the cutoff , sco , and 
 0 otherwise . like the density cutoff , gco , the values of eco and sco are fitted parameters and hence are free to go to zero 
 if imposing these cutoffs does not improve the accuracy of the scoring 
 function .
the total energy and entropy were computed as eknorm  ek , swnorm + 2ek , wwnorm  2eww ,
the quantity eknorm is the mean interaction energy of the water in voxel k with the protein and all other waters , relative to what its interactions 
 would be in bulk , 2eww , bulknorm , computed in the same convention 
 as the other gist quantities . with the voxels binary 
 masks in place for density , energy , 
 and entropy
, we now define the three candidate scoring functions , 
 one using both the energy and entropy data from gist , the second using 
 only the energy data , and the third using only the entropy data:1 here , eaff and saff are additional 
 fitted parameters , which specify the affinity 
 increments provided by voxels surpassing the energy and entropy cutoffs , 
 respectively , and also meeting the criteria dk = gk = 1 .
note that in this study , each of these scoring functions 
 was trained separately and has its own fitted values of gco and c , as well as eco and eaff and/or sco and saff .
we adjusted the scoring functions described above 
 to fit the measured relative binding free energies , gexpt , of 28 different congeneric pairs of fxa 
 inhibitors ( see below ) .
separate training and test sets were used , 
 in order to avoid overfitting of the parameters .
thus , we used a random 
 number generator to split the 28 pairs into two arbitrarily selected 
 sets of 14 apiece .
ten such random splits were carried out , creating 
 10 distinct training and test sets .
parameters were optimized for 
 each training set and then tested on the corresponding test set .
we 
 report means and standard deviations over these 10 splits for the 
 resulting fitted parameters and accuracy metrics .
we further assess 
 the significance of these results by comparing them with results obtained 
 after a shuffling operation , which used the gsl_permutation function 
 in the gnu scientific library , to randomly exchange the entropy and 
 enthalpy values among pairs of voxels . for each training set , 
 the parameters were adjusted as follows , using ges as an example .
we scanned values 
 of eco and sco from 0 to 4 kcal / mol and values of gco from 0 to 4 ( in units of neat water density ) , each in increments 
 of 0.1 .
this scan yields 41  41  41 = 68  921 combinations 
 of the three cutoff values . for each combination ,
linear regression 
 was then used to obtain values of eaff and saff that provide the highest correlation 
 coefficient ( r ) of the relative scores 
 for the congeneric pairs , ges , to the corresponding experimental values , gexpt , for the training set .
the optimized values 
 of the five fitted parameters were then used to compute ges for the congeneric pairs in the training 
 set , and the reliability of the scoring function was evaluated based 
 on the resulting value of r for the test 
 set .
analogous procedures were used for the other two scoring functions , 
 ge and gs .
these require scanning only 
 41  41 = 1681 cutoff combinations and yield only one of the 
 two affinity parameters , eaff or saff , rather than both as for ges .
hydration sites in the fxa binding site were 
 defined and analyzed thermodynamically based on the same md simulation 
 used for the gist calculations , using the first 10 ns ( 10  000 
 frames ) , in accord with the common practice of using approximately 
 210 ns simulations for hsa calculations .
we first collected 
 all instances , in these 1000 frames , of water molecules within 5  
 of any heavy atom of any bound ligand ( see below ) .
for each water 
 molecule in this set , we counted the number of neighboring waters 
 from the same set , using the criterion of an oxygen 
, a water molecule can 
 count as its own neighbor , if two instances of it in different frames 
 meet the distance criterion .
the location of the first hydration site 
 was then set to the coordinates of the water oxygen with the most 
 neighbors .
this water molecule and all of its neighbors were then 
 removed from consideration as potential hydration sites , and the location 
 of the next hydration site was set to the coordinates of the remaining 
 water oxygen with the most neighbors , based on the initial counts . 

this removal process was iterated until the number of neighbors of 
 all remaining waters was less than twice that expected for a 1000 
 frame simulation of bulk water ( i.e. , < 280 from 1000 frames ) . 

each hydration site then was associated with all water instances , 
 from the full 10  000 md frames ( above ) , whose oxygens lay within 
 1  of the site .
each hydration site i was associated with a mean energy ei and a one - body entropy si .
the energy of a water molecule in a given hydration 
 site was calculated as half the difference between the total energy 
 of the water 
a script invoking the program amber , with settings matched to those of the md simulation , was used to 
 compute these energies .
the mean energy of the hydration site then 
 is the average of these energies for all water molecules that populate 
 the site , minus the average energy of a water molecule in neat water 
 from matched calculations .
the water entropy si associated with hydration site i is the sum of its one - body translational and orientational entropies , si , swtrans 
k g(r ) ln g(r ) dr and si , sworient  ( ( kni)/ ) g( ) ln g( ) d , where r is position in 
 the protein frame of reference , k is boltzmann s 
 constant ,  is bulk water density , g(r ) is the local water density referenced to bulk , ni is the number of water molecules 
 associated with hydration site i ,   
 8 , vhs indicates an 
 integral restricted to the spherical hydration site , and  
 defines orientational coordinates in the protein frame of reference . 

the translational entropy was computed by the histogram method , where 
 spherical coordinates r ,  ,  centered 
 on the hydration site were divided into uniformly spaced bins in r , cos  ,  to generate 512 three - dimensional 
 bins of equal volume , with r  [ 0 , 1 ]  , 
   [ 0 ,  ] , and   [ 0 , 2 ] . 

the orientational entropy associated with a hydration site was computed 
 via the same nearest neighbor method used by gist for individual voxels .
we used 
 the hydration sites described above as the basis for three cutoff - based 
 scoring functions , whose functional forms build on prior work . like the three gist scoring functions ( above ) , 

the free hsa - based scoring functions are based on , respectively , both 
 energy and entropy , energy alone , and entropy alone:2 here , the sums 
 range over hydration 
 sites i ; ei and si equal 
 1 if ei and tsi are greater than cutoff values eco and sco , respectively , 
 and 0 otherwise ; di is 
 a displacement function , defined below , which accounts for the overlap 
 of the ligand , in a given pose , with hydration site i ; eaff and saff are fitted constants ; and c is a constant offset . 

note that the hsa scoring parameters in eq 2 are set independently of the gist scoring parameters in eq 1 , despite the use of some equivalent symbols .
one is identical 
 to that used in the previous paper , while 
 the second , as discussed below , applies a physically motivated cap 
 to this quantity:3 here , (x ) , 
 the heaviside step function , 
 equals 0 if x  0 and 1 otherwise ; rco is a distance cutoff ; rij is the distance between hydration site i and atom j of the ligand being scored ; 
 and the sum runs over all atoms belonging to the ligand being scored , i. for dinocap , which is modeled on the prior 
 hydration site scoring function , each 
 ligand atom within rco of the hydration 
 site makes a contribution that scales between 0 and 1 as it approaches 
 the center of the site , so the displacement function accounts for 
 the degree to which a ligand displaces the water in the site .
however , 
 this approach is nonphysical in the sense that , if rij < rco for multiple ligand atoms j and a single site i , then the displacement of solvent from the site may be 
 multiply counted .
that is , for the energy / entropy scoring function , 
 a site might contribute more than eaff + saff ; for the energy - only scoring 
 function , it might contribute more than eaff ; and for the entropy - only scoring function , it might contribute 
 more than saff . indeed , in the present 
 study , we found values of dinocap up to 3.2 , implying 
 that a single site might contribute over three times .
we therefore 
 also considered the alternative displacement function , dicap , which is capped at a value of 1 , so that no hydration 
 site may contribute more than one - fold to a ligand s score , 
 in accordance with the fact that the site can not be displaced more 
 than once .
the parameters of the hsa - based scoring functions 
 were adjusted in the same manner as those of the gist - based scoring 
 functions , except that rco took the place 
 of gco .
thus , values of rco were scanned from 2 to 3  , in steps of 0.1  
 and , as for gist , values of eco and sco were scanned from 0 to 4 kcal / mol in increments 
 of 0.1 kcal / mol .
( note that no hydration sites had energies greater 
 than 3.7 kcal / mol . )
the sums in eq 2 were evaluated 
 for each ligand and with each combination of rco , eco , and sco . for each set of cutoffs scanned ,
values of eaff and saff were 
 obtained by linear regression against the differences in measured 
 binding free energies for a training set of congeneric ligand pairs , 
 for each scoring function in eq 2 .
the parameters 
 that yielded the highest correlation coefficients were chosen and 
 were tested for their ability to reproduce the difference in binding 
 affinities for the congeneric pairs in the test set .
this procedure 
 was applied to the same 10 training and test sets used for the gist 
 scoring function . as done for gist , we assessed the significance of 
 the hsa results by comparing them with results obtained after a shuffling 
 operation , which randomly exchanges the entropy and enthalpy values 
 across pairs of hydration sites . both gist and
the hsa methods take as input a 
 boltzmann sample of water configurations for a given configuration 
 of the protein . here , we used molecular dynamics ( md ) simulations 
 to generate this sample from the canonical ensemble ( nvt ) , as follows . 

we used the structure of fxa from protein data bank entry 1fjs , as previously done , and our assignment of protonation states was also consistent 
 with this prior study .
we removed the ligand from the binding site 
 and used tleap and other amber tools to 
 assign protein parameters from the amber99sb force field and solvate the protein with 8557 tip3p water 
 molecules .
the simulations used a periodic 
 box with dimensions 66.5   72.2   60.9  , 
 which afforded at least 10  between any protein atom and the 
 edge of the periodic box .
four disulfide bonds were set up to join 
 cysteine pairs 7/12 , 27/43 , 156/170 , and 181/209 , and two ions observed 
 in the crystal structure ( ca and cl ) were restrained to their original positions .
the resulting simulation 
 system had 29  338 atoms , comprising the protein , the two ions , 
 and the water molecules .
energy minimization , followed by md 
 simulation , was carried out with the amber 12 software using pmemd.cuda on a single gpu .
first , the energy of the system 
 was minimized in two rounds ; both used 1500 steps of the steepest 
 descents algorithm followed by the conjugate gradient method for a 
 maximum of 2000 steps . in the first round , all protein atoms were 
 harmonically restrained to their initial positions with a force constant 
 of 100 kcal / mol /  . in the second round , the system 
 was further relaxed keeping only non - hydrogen protein atoms restrained , 
 with the same force constant .
the energy minimized system was then 
 heated with a series of 20 ps constant - volume and -temperature md 
 simulations with the first simulation at 50 k and the temperature 
 incremented by 50 k every 20 ps until 300 k was reached .
the system 
 was then equilibrated for 10 ns at 300 k at a constant pressure of 
 1 atm . at the final volume
, the system was then equilibrated for an 
 additional 5 ns at constant volume .
the final md production run of 
 100 ns was at constant number of particles , volume , and temperature 
 ( nvt ) , and system configurations were stored every 1 ps , for a total 
 of 100  000 stored configurations . during all md simulations , 

all protein atoms were harmonically restrained to their positions 
 following the energy minimization step , with a force constant of 100 
 kcal / mol /  .
the shake algorithm was 
 used to constrain the lengths of all bonds involving hydrogen atoms . 

particle mesh ewald was implemented to account 
 for long - range electrostatic interactions , and the leapfrog algorithm 
 was used to propagate the trajectory . for the constant pressure simulations , 
 isotropic position scaling was implemented with a pressure relaxation 
 time of 0.5 ps . the main gist and hsa solvation maps
however , in order to study the convergence 
 properties of gist , we performed two additional 20 ns nvt production 
 runs storing configurations every 0.05 ps ; one was begun identically 
 with the 100 ns run , and the other was begun with the last md configuration 
 of the 100 ns run .
we trained and tested the scoring functions with a set of 28 congeneric 
 ligand pairs ( see supporting information ) , where the members of each pair differ only by small , localized 
 chemical changes in rigid moieties , leading to differential displacement 
 of solvent .
as previously discussed , this 
 approach minimizes the contributions of free energy terms other than 
 hydration and thus allows a focus on the quantity of central interest 
 in this study .
the 28 pairs used here are a subset of 31 drawn from 
 several experimental series for use in a previous computational study : we eliminated three pairs ( matter:25/matter:28 ; matter:28/2bmg : i1h ; 
 mueller:3/mueller:2 ) , because we were not confident of the conformation 
 of at least one member of the pair , and hence of the location of the 
 displaced solvent .
in particular , mueller:3 differs from mueller:2 
 by a phenyl group whose orientation is not clear , because it is attached 
 by a rotatable bond ; and for the other two pairs , an aromatic ring 
 changes to a nonaromatic ring , whose conformation is uncertain . it 

is worth noting that , for this set of ligands , there is essentially 
 no correlation between binding free energy and molecular weight ( r = 0.12 ) .
ligand poses were drawn from 
 available cocrystal structures or generated from a cocrystal structure 
 of a closely related ligand by a small chemical adjustment . in all 
 cases , the cocrystal structure was aligned with the simulated protein 
 ( above ) to generate an initial pose in the binding site for which 
 the hydration structure was computed .
final poses for modeling solvent 
 displacement were generated by protonating the ligands , then minimizing 
 the initial poses in the simulated protein structure while allowing 
 the ligand and protein hydrogen atoms to relax .
the atomic partial 
 charges for each ligand were obtained from the restrained electrostatic 
 potential ( resp ) method , using quantum - mechanically 
 derived electrostatic potentials with the 6 - 31 g * basis set .
note that these parameters were used only to 
 generate the ligand poses studied with the gist and hsa scoring functions .
the hsa hydration sites are informative 
 about water density but do not capture the level of detail available 
 from gist s grid - based method .
the two approaches are compared 
 in figure 1 , which displays the hsa sites ( blue 
 spheres ) computed for the binding pocket of factor xa , along with 
 gist s gridded representation of water density , contoured at 
 three different levels .
contours of water density at 6 times that 
 of bulk ( g = 6 ) appear as discrete , mostly convex 
 droplets ( yellow contours , left panel ) , although a few of these high - density 
 regions are elongated , rather than round .
every high - density droplet 
 is matched by a spherical hsa site , but there are many hsa sites in 
 the binding pocket that do not enclose one of these high - density droplets . 

contours at 4 times bulk density ( g = 4 ) appear as 
 more and larger droplets and match the hsa sites rather well ( gold 
 contours , middle panel ) . contours at twice bulk density ( g = 2 ) form long , curved strands , which follow the contours of the 
 protein surface ( orange contours , right panel ) , and are not well represented 
 by the hsa sites . these begin to delineate the first hydration shell 
 of the protein .
finally , contours at still lower density ( e.g. , g = 1.5 ) include parts of the second hydration shell ( not 
 shown ) .
overall , the hsa representation captures the droplet - like 
 distribution of the highest water densities but does not distinguish 
 between high and medium density regions and does not capture the complex 
 distributions of water density that become apparent at densities roughly 
 twice that of bulk .
this observation has practical relevance , because 
 the scoring functions developed here include regions where water density 
 is of this order , as detailed in the next subsection .
comparison of the gist 
 and hsa representations of water density 
 in the factor xa binding pocket .
hsa sites ( blue spheres ) are the 
 same in all three panels . from left to right , the gist contour levels 
 are at g = 6 , 4 , and 2 .
the gist water densities 
 are based on the occupancy of grid voxels by water oxygens , and the 
 boundaries of the grid box may be discerned in the right - hand panel . 

it is also worth remarking that 
 a hydration site should not be 
 directly equated with a single bound water , as the sites occupancies 
 are typically well below one .
thus , for the hsa sites within 5  
 of the bound ligands , we find a mean occupancy of 0.58 waters molecules , 
 with a standard deviation of 0.22 .
these values are similar to those 
 reported previously for factor xa in an early implementation of watermap : the occupancies reported in table 1 of the prior 
 report correspond to a mean of 0.51 water molecules per hydration 
 site , with a standard deviation of 0.22 . a gist - based 
 scoring function based on both the local energy and one - body entropy 
 of displaced water yields good correlations with the measured binding 
 affinity differences for the 28 congeneric ligand pairs ( table 1 , top row ) .
the mean r value is 0.94 across 10 training sets drawn at random from the full 
 set of 28 pairs , and a high correlation ( r = 0.84 ) is preserved when the optimized parameters are applied to 
 the respective test sets .
the best results are obtained when voxels 
 are excluded if their water density is less than about twice the bulk 
 density ( gco = 1.9 ) .
each voxel which 
 furthermore meets the energy cutoff ( about 0.6 kcal / mol / water above 
 bulk ) contributes about 0.2 kcal / mol of affinity .
the entropy 
 terms are somewhat puzzling : the entropy cutoff of 3.3 kcal / mol / water 
 is much higher than the energy cutoff of 0.6 kcal / mol / water , and the 
 positive value of saff seemingly indicates 
 that displacing low entropy water disfavors binding , rather than favoring 
 it , as might be anticipated .
these results suggest that the energy 
 term may be more meaningful ; this topic is further discussed below .
the quality 
 of each scoring function 
 is reported in terms of coefficients of determination ( r ) .
the test - set results ( test r ) are the most meaningful ; the training - set results ( train r ) are included for completeness .
all quantities 
 are averages across the 10 splits , with associated standard deviations . 

results are presented for the actual hydration data and for hydration 
 data shuffled among voxels ( gist ) or hydration sites ( hsa ) , and hsa 
 results are presented for capped and uncapped displacement functions 
 ( see methods ) .
scoring functions were constructed 
 based on both energy and entropy data ( e / s ) , energy data only ( e ) , and entropy data 
 only ( s ) ; the parameters are defined in the methods section .
it is of interest to visualize the energy and entropy 
 scoring regions ; 
 i.e. , those voxels which meet either both the density and energy cutoffs 
 or both the density and entropy cutoffs , respectively .
as shown in 
 figure 2 ( left ) , the energy scoring region 
 tends to localize at extended regions of the nonpolar surface .
these 
 are places where water molecules lose more energy due to desolvation 
 by the protein than they gain by making new interactions with the 
 protein .
( the projection from three to two dimensions makes some scoring 
 regions appear deceptively close to polar protein atoms . )
the analogous 
 visualization for the energy - only scoring function , which is discussed 
 below , is very similar ( data not shown ) .
the locations of the entropy 
 scoring regions ( figure 2 , middle ) are more 
 complicated .
some ( white arrows ) lie at the surface of polar atoms ; 
 others ( pink arrows ) lie at hydrophobic locations .
the frequent localization 
 of entropy scoring regions at polar surfaces , and the unexpectedly 
 positive value of saff , likely reflects 
 the fact that polar surfaces can tightly bind waters , leading to unfavorable 
 entropies but favorable energies ( figure 2 , 
 right ) typical of traditional entropy  enthalpy compensation . 

displacement of water from such regions may be net unfavorable , and 
 this might help account for the positive value of saff . on the other hand , the displacement of entropically 
 disfavored water from subregions where energy is not particularly 
 favorable should favor ligand binding .
thus , the energetically mixed 
 nature of the entropic scoring regions further suggests that it may 
 not give a clear signal in the overall scoring function .
gist contours 
 in binding site of fxa , in molecular surface representation . 
 left : energy scoring region , which meets both density and energy cutoff 
 criteria for the combined energy 
middle : entropy scoring region , which meets both density and entropy 
 cutoff criteria for the combined energy  entropy scoring function . 

right : normalized total water energy , contoured at 2.6 kcal / mol / water . 

graphic generated with vmd . in order to study the significance of the energy and entropy 
 terms 
 in more detail , we also considered a scoring function based only on 
 density and energy , and another based only on density and entropy . 

as shown in table 1 ( second and third rows ) , 
 the energy - only scoring function performs just as well ( r = 0.85 for the test sets ) as the original one based 
 on both energy and entropy .
in addition , the fitted parameters are 
 similar in magnitude and sign to those of the original energy term . 

in contrast , the entropy - only scoring function yielded a poor correlation 
 with experimental data ( r = 0.35 for 
 the test sets ) , and the sign of saff is 
 reversed relative to that in the original scoring function .
thus , 
 the hydration energy alone carries all of the predictive power of 
 the gist - based scoring function , at least for fxa with this ansatz . 

this result is consistent with the analysis of the combined energy / entropy 
 scoring function , above . as a further check of the statistical 
 significance of the present 
 results
, we shuffled the gist data among voxels and then refitted 
 all three gist - based scoring functions using the shuffled data .
although correlations as high as r = 
 0.82 are obtained for the training sets , the test - set results are 
 all poor ( r  0.4  0.2 ) . 

this result supports the significance of the high correlations obtained 
 with the true ( unshuffled ) data and indicates that the low correlations 
 observed for the entropy - only scoring function should be viewed as 
 statistically insignificant . finally , we examined the amount 
 of md simulation data required 
 to generate the high correlations observed above .
first , we reanalyzed 
 the original set of md frames , which had been saved at 1 ps intervals . 
 for increasing numbers of frames from this set , we reran the 10 training 
 and test calculations and computed the mean and standard deviation 
 of the resulting 10 values of r.
as shown 
 in figure 3 ( left ) , the value of r and its standard deviation ( error bars ) appear to plateau 
 at about 60 ns for the combined energy  entropy scoring function 
 ( top left ) .
interestingly , the plateau starts much earlier , at about 
 30 ns , for the energy - alone scoring function ( lower left ) .
it appears 
 that the more slowly convergent entropy term delays convergence of 
 the energy 
entropy scoring function , so that leaving out the 
 entropy term in the energy - only scoring function speeds convergence . 

this result is , again , consistent with the irrelevance of the entropy 
 term in the gist scoring functions .
we then asked whether shorter 
 md simulations might give better convergence if frames were saved 
 at shorter time intervals .
first , we reprocessed the first 20 ns of 
 the same simulation , now processing frames saved at 0.05 ps . as shown 
 in figure 3 ( top middle ) , the combined energy / entropy 
 scoring function now converges somewhat more quickly , especially as 
 measured by the reduction in the standard deviation of r across the 10 train / test calculations .
the improvement 
 is more marked for the energy - only scoring function , as well - converged 
 results are now available after only 10 ns of simulation time , although 
 the standard deviation of r remains slightly 
 higher than that from the 100 ns simulation .
we then extended the 
 100 ns simulation by 20 ns , saving frames at 0.05 ps intervals , and 
 examined convergence over this short simulation .
the results are further 
 improved , with good convergence and tight error bars achieved within 
 about 5 ns for the energy - only scoring function , and 10 ns for the 
 combined energy / entropy one .
the improvement in these results , relative 
 to those from the early 20 ns segment , suggests that the water structure 
 continued to equilibrate somewhat during the 100 ns run , so it would 
 have been appropriate to use a somewhat longer equilibration period 
 in the md protocol . in summary , at least for fxa , a simulation of 
 10 ns or less suffices to gain all the benefit of these gist scoring 
 functions .
convergence of r values for gist scoring 
 functions , as a function of simulation duration .
bottom row : energy - only scoring function . 
 left : 100 ns simulation , frames saved at 1 ps intervals .
middle : the 
 first 20 ns of the same 100 ns simulation , frames saved every 0.05 
 ps .
right : 20 ns simulation initiated from the last frame of the 100 
 ns simulation , frames saved every 0.05 ps . as discussed 
 in the methods section , a prior hsa scoring 
 function was constructed in such a way that the thermodynamic contribution 
 of each hydration could be counted multiple times , if more than one 
 ligand atom lay within a cutoff distance .
this function gave good results but is arguably nonphysical , because 
 once the water in a hydration site has been displaced , it can not be 
 displaced again .
here , we present results for a set of similarly constructed 
 hsa - based scoring functions in which , as before , a hydration site 
 can be counted multiple times ; as well as a second set , in which the 
 contribution of each hydration site is capped , so that it can only 
 contribute once .
for both the uncapped and capped models , we examine 
 scoring functions based on energy and entropy , energy only , and entropy 
 only , as also done for gist . with the original uncapped approach , 
 where a hydration site can contribute multiple times
, the combined 
 energy / entropy scoring function provides high correlation coefficients 
 for both the training ( r = 0.88 ) and 
 test ( r = 0.80 ) sets ( table 1 , fourth row , marked e / s ) .
these values of r are slightly 
 lower than those for gist , but the difference is not statistically 
 significant .
interestingly , the fitted value of eaff ( 1.70 kcal / mol ) is much greater than that 
 of saff , which in fact is assigned a positive 
 sign ( 0.09 kcal / mol ) .
in addition , the training procedure puts much 
 sharper constraints on the energy terms than the entropy terms , as 
 indicated by the fact that the standard deviations of eo and eaff ( 0.08 and 0.31 
 kcal / mol , respectively ) are lower than those for so and saff ( 0.50 and 0.60 
 kcal / mol ) . thus , the training results suggest that the entropy term 
 is of lower importance than the energy term , much as observed in the 
 case of gist . accordingly , a scoring function based entirely on energy 
 ( table 1 , row 5 , marked e ) 
 performs as well as the one with both energy and entropy , and furthermore 
 yields values of rco , eco , and eaff very similar 
 to those of the combined energy  entropy scoring function . on 
 the other hand , an entropy - only scoring function ( table 1 , row 6 , marked s ) also performs fairly well , with a value 
 of r = 0.66 for the test sets , and the 
 fitted value of saff , 1.52 kcal / mol , 
 is essentially the same as the fitted value of eaff for the energy - only scoring function , 1.50 kcal / mol . 

thus , the entropy - only scoring function appears to largely replicate 
 the energy - only scoring function , with some drop in the correlation 
 with experimental data .
these results are similar to those for gist , 
 as in both cases , the energy - only scoring function performed as well 
 as the energy / entropy one , while the entropy - only one is worse . however , 

the decline in performance on going to entropy - only is much greater 
 for gist than for hsa . in the second variant of the hsa - based 
 scoring functions ,
no hydration 
 site can contribute more than one - fold to the difference between two 
 ligands affinities . imposing this physically reasonable
cap 
 on the contribution of each site reduced the experimental correlation 
 of the scoring function with experimental data to 0.66 for the test 
 sets with the combined energy 
entropy scoring function but 
 had essentially no effect on the correlations for the energy - only 
 and entropy - only test - set results .
thus , the results remain consistent 
 with a conclusion that the energy term alone is enough to gain all 
 the benefit of the scoring functions .
the only other major change , 
 relative to the uncapped version of the scoring function , is that 
 the fitted values of eaff and saff changed from about 1.5 kcal / mol 
 to about 2.3 kcal / mol , except that saff for the combined energy / entropy scoring function remained 
 small .
these increases presumably have the effect of compensating 
 for the reduced values of the displacement function , i.e. , for the 
 fact that dicap  dinocap in eqs 2 and 3 . in order to further test the statistical significance of the hsa 
 scoring results , we shuffled the energies and entropies among hydration 
 sites and then refitted all six hsa - based scoring functions with the 
 shuffled data .
as summarized in the last six rows of table 1 , the uncapped scoring functions now yield poor 
 correlations with experimental data ( r  0.4 ) , just as observed for gist , and the correlations for 
 the capped scoring functions fall even lower .
these results support 
 the significance of the correlations obtained with the actual ( unshuffled ) 
 data .
the fact that the shuffled results are worse for the capped 
 hsa scoring functions suggests , but does not prove , that applying 
 the physically reasonable cap might reduce spurious correlations .
it is of interest to examine the locations of hsa scoring sites 
 relative to the gist scoring regions .
we focus here on the energy - only 
 scoring functions , since it is not clear that including entropy adds 
 useful information , at least for fxa within the present functional 
 form . as shown in figure 4 ( left ) , the hsa 
 sites do not capture the complex shapes of the gist scoring regions , 
 and substantial parts of the gist scoring region are entirely missed 
 by the hsa scoring sites . these regions
the spatial relationship of the 
 scoring sites to the ligands studied here may be appreciated from 
 figure 4 ( right ) , which shows the van der waals 
 surface of a representative ligand .
comparison of gist and hsa scoring locations , 
 in the context of 
 the fxa binding site ( gray ) . left : scoring sites from the capped hsa 
 energy - only scoring function ( blue spheres ) and energy scoring regions 
 from the gist energy - only scoring function ( orange contours ) .
gist 
 results are present only within the limits of the gist grid ( orange 
 lines ) .
right : a representative ligand , pdb het i d 4qc , in van der 
 waals representation ( pink ) .
ligand 
 modeling is that the detailed representation of water structure and 
 thermodynamics it affords works at least as well in a simple scoring 
 function as the prior and present site - based 
 hsa implementations .
in addition , the gist results converge well within 
 510 ns of md simulation time , depending upon whether one uses 
 the energy - only model or the energy / entropy model .
these simulation 
 times are commensurate with those normally used for the hsa approach ; 
 see the methods .
additional 
 considerations include the fact that site - based approaches may require 
 less data storage and that they paint a simple picture of water structure 
 in a protein binding site . on the other hand ,
the gist grid files 
 produced are still small ( < 1 mb ) , and we anticipate that the more 
 detailed rendering of hydration structure and thermodynamics afforded 
 by the grid approach will be useful for insight and prediction .
note , 
 in particular , that there are hsa sites which are only partly occupied 
 by the more refined gist scoring regions , as well as gist scoring 
 regions that are not identified by the hsa model , as evident from 
 figure 4 .
more generally , the fact that gist 
 is a more direct representation of the common underlying inhomogeneous 
 solvation theory facilitates the interpretation of its results and 
 provides clear pathways to future enhancements , such as the incorporation 
 of higher - order correlations , as touched on below .
a fast implementation 
 of gist , based upon focused grand canonical monte carlo sampling , may be of particular interest in the near term 
 for high - throughput applications .
another novel and striking 
 result of the present study is that 
 both the gist and hsa models provide clear signals that ligands can 
 gain affinity by displacing energetically unfavorable binding site 
 water , whereas the displacement of entropically unfavorable water 
 seems to play a negligible role .
the energetically unfavorable water 
 highlighted by this study localizes at hydrophobic patches of the 
 protein surface , perhaps especially in concave regions where water 
 molecules are expected to lose hydrogen bonds .
the concept that energy may outweigh entropy in cases of 
 strong hydrophobic binding has been raised before , in both experimental and computational contexts . nonetheless , water entropy is typically thought to play 
 a central role in hydrophobic binding . here ,
interestingly , neither 
 the gist or hsa models made a compelling case that the displacement 
 of entropically unfavorable water consistently enhances affinity . 

we conjecture that the lack of a clear correlation of water entropy 
 with affinity may reflect the fact that low water entropy often results 
 from energetically favorable water 
protein interactions , so 
 that water may actually be quite stable in many locations where its 
 entropy is low .
this view is consistent with the experimental observation 
 that the entropy of hydration of small ions is strongly negative , although the free energy is also strongly negative . 

it may be possible in the future to devise a more sophisticated scoring 
 function that would account for the enthalpy  entropy compensation 
 between stabilizing energy and destabilizing entropy and focus on 
 regions where this compensation breaks down , such as in the binding 
 cavity of the synthetic host molecule cucurbituril .
it is also worth mentioning that different protein 
 binding sites affect water differently , so a different result might 
 be obtained for a different protein .
finally , it may be that capturing 
 the entropic aspect of the hydrophobic effect requires accounting 
 for water  water correlations , which are absent from the one - body 
 entropy considered here .
if so , the entropy term may become more important 
 once pairwise correlations have been incorporated into gist s 
 entropy calculations .
our observation that the displacement 
 of high energy water plays 
 a greater role in ligand scoring than displacement of low entropy 
 water appears to contrast with a prior hsa - based study of the same 
 system , where the fitted scoring function placed approximately equal 
 weight on water energy and entropy .
however , 
 the range of fitted values for saff in 
 our hsa models nearly spans the value of 0.66 kcal / mol for the corresponding 
 parameter in the prior study , srwd . in 
 addition
, the prior study did not examine the uncertainty in its fitted 
 energy and entropy parameters or evaluate a scoring function based 
 purely on water energy .
therefore , the results of these two studies 
 should not be regarded as inconsistent . in summary ,
the grid - based 
 gist method of extracting information 
 about hydration thermodynamics from md simulations with explicit water 
 has provided encouraging results in its first application to protein  ligand 
 binding .
it thus appears to hold significant promise as a broadly 
 applicable method of understanding the role of binding site water 
 in protein  ligand binding , and as a tool to improve the accuracy 
 of methods for discovering high affinity targeted ligands .
it is anticipated 
 that the detailed representation of water distributions and thermodynamics 
 which gist affords will make it particularly informative .
we are currently 
 working to develop such applications and to release an open - source 
 implementation of gist within the ambertools software package for others to study and use .
 water molecules in the active site 
 of an enzyme occupy a complex , 
 heterogeneous environment , and the thermodynamic properties of active - site 
 water are functions of position . as a consequence , it is thought that 
 an enzyme inhibitor can gain affinity by extending into a region occupied 
 by unfavorable water or lose affinity by displacing water from a region 
 where it was relatively stable . 
 recent advances in the characterization 
 of binding - site water , based on the analysis of molecular simulations 
 with explicit water molecules , have focused largely on simplified 
 representations of water as occupying well - defined hydration sites . 
 
 our grid - based treatment of hydration , gist , offers a more complete 
 picture of the complex distributions of water properties , but it has 
 not yet been applied to proteins . 
 this first application of gist to 
 protein  ligand modeling , for the case of coagulation factor 
 xa , shows that ligand scoring functions based on gist perform at least 
 as well as scoring functions based on a hydration - site approach ( hsa ) , 
 when applied to exactly the same simulation data . 
 interestingly , the 
 displacement of energetically unfavorable water emerges as the dominant 
 factor in the fitted scoring functions , for both gist and hsa methods , 
 while water entropy plays a secondary role , at least in the present 
 context . 
Introduction
Methods
Results
Discussion