  a high - throughput virtual screening pipeline has been extended from single energetically minimized structure molecular mechanics / generalized born surface area ( mm / gbsa ) rescoring to ensemble - average mm / gbsa rescoring .
the correlation coefficient ( r2 ) of calculated and experimental binding free energies for a series of antithrombin ligands has been improved from 0.36 to 0.69 when switching from the single - structure mm / gbsa rescoring to ensemble - average one .
the electrostatic interactions in both solute and solvent are identified to play an important role in determining the binding free energy after the decomposition of the calculated binding free energy .
the increasing negative charge of the compounds provides a more favorable electrostatic energy change but creates a higher penalty for the solvation free energy .
such a penalty is compensated by the electrostatic energy change , which results in a better binding affinity .
a highly hydrophobic ligand is determined by the docking program to be a non - specific binder .
  our results have demonstrated that it is very important to keep a few top poses for rescoring , if the binding is non - specific or the binding mode is not well determined by the docking calculation .
high - throughput virtual screening is one of the most commonly used techniques in the drug discovery procedure .
a high throughput virtual screening pipeline has been built to leverage the high performance computing ( hpc ) resource at lawrence livermore national laboratory for in - silico screening of large virtual compound databases .
the pipeline consists of five modules : receptor / target preparation , ligand preparation , vinalc docking calculation , single - structure mm / gbsa rescoring , and ensemble - average mm / gbsa rescoring .
the pipeline uses a down - selecting scheme to filter large numbers of compounds through the docking , single - structure rescoring , and ensemble - average rescoring .
one of the notable features of this pipeline is an automated receptor preparation scheme with unsupervised binding site identification , which enables automatically running the whole pipeline with little or no human intervention .
perez - sanchez and co - workers have developed a similar approach to improve drug discovery using massively parallel gpu hardware instead of supercomputers .
their gpu - based program , bindsurf , takes advantage of massively parallel and high arithmetic intensity of gpus to speed - up the calculation in low cost desktop machine . in molecular docking ,
the scoring functions simplify the calculation with many approximations to achieve high throughput screening . in the previous study
, we have demonstrated the parallel vinalc docking program can screen 1 million compounds in less than 1.5 hours by using 15k cpus .
docking programs sacrifice accuracy for speed to make the assessments of a large number of compounds feasible . even though the ligand binding conformation can be identified in most cases by using various existing docking programs , there are no universal scoring functions for all types of molecules and protein families .
rescoring steps are necessary for the accurate prediction of binding energies . by using approximated scoring functions , molecular docking usually fails to estimate binding energies when compared to experimental values .
the mm / gbsa method is selected for rescoring because it is the fastest force - field based method that computes the free energy of binding , as compared to the other computational free energy methods , such as free energy perturbation ( fep ) or thermodynamic integration ( ti ) methods .
comparison studies have also shown that mm / gbsa outperforms molecular mechanics / poisson boltzmann surface area mm / pbsa .
the mm / gbsa method has been widely exploited in free energy calculations [ 10 , 11 ] , and mm / gbsa rescoring overly yields better results than docking for the dud - e ( directory of useful decoys , enhanced ) data set [ 1 , 12 ] .
there is evidence showing that molecular dynamics ( md ) simulations are necessary for some systems to identify the correct binding conformations . in this study , we have extended our pipeline from single - structure to ensemble - average mm / gsba rescoring by carrying out long - time md simulations . to validate the new approach ,
1 ) , including heparin and non - polysaccharide scaffold compounds . for the purpose of comparison ,
both single - structure and ensemble - average mm / gsba rescoring are employed in the binding affinity calculations of antithrombin ligands .
we must point out that estimation / calculations of the entropy term are tricky . in most scenarios ,
quite a few researchers dispute the benefits of including the entropy term , which can be a major source of error due to the drawback of the entropy calculation method [ 13 , 14 ] , despite others who advocate its usage .
we choose to neglect the partial entropy effects estimated by the normal mode analysis in our calculations .
antithrombin is a glycoprotein that plays a crucial role in the regulation of blood coagulation by inactivating several enzymes of the coagulation system and , thus , is an important drug target for the anticoagulant treatment .
-antithrombin is the dominant form of antithrombin and consists of 432 amino acids with 4 glycosylation sites , where an oligosaccharide occupies each glycosylation site .
heparin is the first compound that was identified and used as an anticoagulant and antithrombotic agent .
1 , ntp ) that binds and activates the antithrombin . this binding localizes the function of antithrombin to inhibition of serine proteases in the coagulation cascade in the bloodstream , which allows coagulant activity in damaged tissue outside the vascular system . due to increasing interests in clinical applications ,
verli and co - workers performed molecular dynamics simulations to study the induced - fit mechanism of the antithrombin - heparin interaction and effects of glycosylation on heparin binding [ 19 , 20 ] .
they also confirmed an intermediate state between the native and activated forms of antithrombin . because of the weak surface complementarity and the high charge density of the sulfated sugar chain
, the docking of heparin to its protein partners presents a challenging task for computational docking .
navarro - fernandez and colleagues screened a large database in silico by using flexscreen docking program   and identified a new , non - polysaccharide scaffold able to interact with the heparin binding domain of antithrombin .
they predicted d - myo - inositol 3,4,5,6-tetrakisphosphate ( fig . 1 , l1c4 ) to strongly interact with antithrombin , which was confirmed by experimental binding affinity study .
here , we carried out molecular dynamics ( md ) simulations of ligand recognition upon binding antithrombin and calculate the ligand binding affinity . according to the characteristics of the compounds , we divide the study into two categories
the first category is devoted to the ligands that specifically binding in the active sites of the antithrombin .
compound l1c1 , l1c2 , l1c3 , l1c4 , l1c5 , l1c6 , and ntp belong to this category .
the binding affinity of the top pose of each ligand - complex is computed by using the ensemble - average mm / gbsa rescoring method , as a complementary study to the previous docking work [ 21 , 23 ] .
the advantage of long - time md simulations is that more configurational space can be explored and dynamical properties of systems can be revealed .
as shown in the following section , penfluridol is a non - specific binder as determined by the docking program .
we selected top 10 poses and three additional poses for ensemble - average mm / gbsa rescoring .
the initial structure for md simulation of the antithrombin complex with compound ntp is obtained from the pdb bank ( pdb i d : 1azx ) .
the initial structures of the antithrombin complex with compound l1c1 , l1c2 , l1c3 , l1c4 , l1c5 , and l1c6 ( fig . 1 ) are obtained from the flexscreen docking program .
the initial structures of the antithrombin complex with penfluridol are the top 10 poses from the docking calculation and three additional poses , where they are close to the active site of antithrombin .
the mm / gbsa calculations are applied to these initial structures by using our in - house developed pipeline [ 1 , 2 ] and amber molecular simulation package .
ligands use the amber gaff forcefield   as determined by the antechamber program   in the amber package .
the fourth module of the pipeline is employed for the single - structure mm / gbsa calculation , where the receptor - ligand complexes are energetically minimized by the mm / gbsa method as implemented in the sander program of the amber package .
the atomic radii developed by onufriev and coworkers ( amber input parameter igb=5 ) are chosen for all gb calculations .
for the ensemble - average mm / gbsa rescoring , a md trajectory is obtained for each ligand - receptor complex .
the snapshots of complex are extracted from the md trajectory at every 10 ps , which yields 10,000 snapshots .
the various free energy terms of complexes , receptors , and ligands derived from snapshots are calculated by the mm / gbsa methods .
the final binding free energy is an average of the binding free energies of snapshot .
the md simulation for each ligand - receptor complex starts with the energetically minimized structure from single - structure mm / gbsa rescoring .
the systems are heated from 0 k to room temperature , 300 k. the md simulations with a time step of 2 fs for the integration of the equations of motion are carried out at room temperature .
binding affinities of antithrombin and its 8 ligands are calculated by post - processing the ensembles of structures extracted from md trajectories using mm / gbsa calculations . in the mm / gbsa calculation , the binding free energy between a receptor and a ligand is calculated using the following equations : the binding free energy ( gbind ) is decomposed into different energy terms . because the structures of complex , receptor , and ligand are extracted from the same trajectory , the internal energy change ( eint ) is canceled .
thus , the gas - phase interaction energy ( egas ) between the receptor and the ligand is the sum of electrostatic ( eele ) and van der waals ( evdw ) interaction energies .
the solvation free energy ( gsol ) is divided into the polar and non - polar energy terms .
the non - polar contribution is calculated based on the solvent - accessible surface area ( gsurf ) .
a value of 80 is used for the solvent dielectric constant , and the solute dielectric constant is set to 1 .
the calculated binding free energy ( gbind ) is the sum of the gas - phase interaction energy and solvation free energy because we neglect the entropy term .
the experimental binding free energy is estimated from the experimental dissociation constant ( kd ) by the equation : where r is the gas constant , and t is temperature . to study the docking and rescoring
, we choose a hydrophobic compound ( penfluridol ) to carry out the docking and the ensemble - average rescoring calculations . to compare docking and rescoring
, we utilize the bindsurf program to explore all possible binding sites on the surface of antithrombin and dock penfluridol to the predicted binding sites by using the flexscreen docking program .
we use bindsurf instead of flexscreen because flexscreen is designed to take into account the structure of the ligand , and the results obtained are similar ( data not shown ) .
the docking results have not shown much difference in docking scores for the top 50 poses ( see results and discussion ) . in the previous study ,
the top 20 docking poses were kept for the mm / gbsa single - structure rescoring .
due to the intensive computational resource requirement in the ensemble - average mm / gsba rescoring calculations , here we keep only the top 10 docking poses for ensemble - average mm / gsba rescoring of one ligand binding study . in the top 10
docking poses , there is only one pose that is docked into the antithrombin active site . because the poses near the active site of antithrombin are our primary interests , we include three additional poses near to the antithrombin active site based on the docking scores .
although these three poses are not in the top 10 list , the docking scores for them are very similar to those in the top 10 list .
the calculated binding free energies of compound l1c1 , l1c2 , l1c3 , l1c4 , l1c5 , l1c6 , and ntp using the ensemble - average mm / gbsa rescoring are shown in table 1 together with their corresponding experimental values .
each calculated binding free energy is averaged from snapshots extracted from 100 ns md trajectories .
as determined experimentally , compound l1c4 is the best binder with a kd value of 0.088 m . as predicted by the mm / gbsa method
, compound l1c4 has the most negative binding free energy ( -308.01 kcal / mol ) , which is in agreement with the experimental results .
the second best binder as predicted by the mm / gbsa calculation is compound ntp with a calculated binding free energy of -279.57 kcal / mol , confirming the experimental ranking relative to compound l1c4 .
compound l1c2 is predicted to have the worst binding free energy of the six ligands , which is also in agreement with its experimental ranking value . in summary ,
the mm / gbsa calculations rank the binding affinities of all six antithrombin ligands in the same exact order as that of experimental binding free energy rankings .
the calculated binding free energies of six antithrombin ligands using the ensemble - average mm / gbsa rescoring have been plotted against the free energies derived from experimental dissociation constants .
the correlation coefficient ( r2 ) is 0.69 , which indicates good correlation between the calculated and experimental values ( fig .
the correlation coefficient calculated by single - structure mm / gbsa is only 0.36 , and compound ntp is predicted to be the best binder , instead of compound l1c4 .
thus , using the ensemble - average mm / gbsa rescoring method significantly improves the accuracy of the prediction over the single - structure mm / gbsa rescoring . as shown in fig .
( 1 ) , compound l1c1 , l1c2 , l1c3 , l1c4 , l1c5 , l1c6 , and ntp contain negatively charged groups , suggesting electrostatic interactions should be a key factor in the binding affinity .
compound ntp has a total charge of -11 , and compound l1c4 has a total charge of -8 .
by decomposing the binding free energy , compound ntp and l1c4 have the largest electrostatic energy changes upon binding in both gas phase ( eele ) and gb solvent ( ggb - ele ) . the energy change upon binding in gas phase is equivalent to the energy change upon binding for the solute .
thus , in other words , compound ntp and l1c4 have the largest electrostatic energy changes upon binding in solute and solvent .
in contrast , compound l1c2 has the smallest electrostatic energy changes in solute and solvent . although compound l1c4 has the least favorable van der waals energy change upon binding , the electrostatic energy change compensates significantly .
for all ligands , the van der waals energy changes ( evdw ) upon binding are less than the electrostatic energy changes ( eele ) by 1 - 2 orders of magnitude .
the contribution of the van der waals energy change has been overpowered by the electrostatic energy change .
non - polar contribution of solvation free energy of compound ntp and l1c3 are more negative than that of the other compounds because the sizes of compound ntp and l1c3 are larger than the other compounds
the non - polar contribution is overwhelmed by the polar contribution of solvation free energy .
thus , the two major factors to determine the binding affinity are the electrostatic energy change and solvation free energy change . the larger the total charge of the compound , the larger the penalty cost is for solvation free energy .
however , the high penalty for large total charge of compound has already been paid by the favorable large changes in electrostatic energy .
although , the electrostatic energy change of compound l1c4 is less than that of compound ntp , compound l1c4 needs less compensation for the solvation free energy .
hydrogen bonding analysis determines the numbers of hydrogen bonds to antithrombin that are persistent at > 20% of the time .
compound ntp has 40 hydrogen bonds to antithrombin , while l1c4 has 25 hydrogen bonds . for compounds l1c1 , l1c2 , l1c3 , l1c5 , and l1c6 , the numbers of hydrogen bonds are 5 , 10 , 12 , 12 , and 12 , respectively .
taking the molecular weight into account and using a similar approach as reynolds ligand efficiency method , compound l1c4 has the highest ligand efficiency .
3a ) . one hydrogen bond ( o6-hh21-nh2 ) is 94.81% persistent , and the other one ( o6-he - ne ) is 89.55% persistent .
the average hydrogen bond distances between the heavy atoms are 2.74  and 2.70  , respectively .
compound l1c4 has strong hydrogen bonds with arg47 , and one of the four phosphate groups from compound l1c4 is locked to the arg47 . according to the hydrogen bonding analysis , compound l1c4 is also hydrogen bonded to arg46 , arg13 , lys114 , lys11 , lys125 , and asn45 , which are key residues to the binding process .
except for the phosphate group locked to arg47 , the other three phosphate groups of compound l1c4 can rotate so that key residues can form hydrogen bonds to different oxygen atoms of phosphate at different times during the md trajectory .
after 8 ns of md simulation , arg13 begins to make hydrogen bonds with the phosphate group of compound l1c4 , which suggests that long - time md simulations are essential to obtaining accurate binding affinities . as shown in fig .
( 3b ) , compound ntp makes hydrogen bonds to antithrombin mainly via its negatively charged sulfate groups .
compound ntp forms high persistent hydrogen bonds with arg13 , arg129 , arg47 , and asn45 ( 70 - 88% ) and forms medium persistent hydrogen bonds with arg132 , lys125 , and thr44 ( 43 - 66% ) .
only relatively weak hydrogen bonds are observed with arg46 , lys114 and lys11 . to compare the docking and rescoring procedures of previous and current studies
the distribution of docking scores shows little energetic difference between the docking poses ( fig .
the top 10 penfluridol docking poses and 3 additional poses near the active site of antithrombin are selected for the mm / gbsa rescoring ( fig .
1 ) compound consists of aromatic rings and alkane groups , which suggests that it is hydrophobic .
thus , the docking sites , as determined by the bindsurf program , are close to hydrophobic surfaces of antithrombin , as shown in the fig .
the active site of antithrombin is highly charged with many positively charged residues , such as arginine and lysine . in the end , most poses have been docked outside of the binding active site of antithrombin .
the ensemble - average mm / gsba rescoring has been performed for the selected 13 docking poses .
the binding free energies for these selected poses range from -28.75 to -59.68 kcal / mol , which are much smaller than those of good antithrombin binders , e.g. compound l1c4 and ntp . as expected , mm / gbsa rescoring ranks the poses differently than the docking method because the docking program uses an empirical scoring function as more like  machine learning 
the scoring function does not contain any exact terms for hydrophobic , vdw , or solvation components . in constrast , mm / gbsa is a physics - based method , which contains explicit terms for hydrophobic , vdw , or solvation components .
the mm / gbsa method is considered to be more accurate at the cost of computer time .
the top docking ranked poses are the lowest ranked poses using mm / gbsa rescoring , indicating the very importance of keeping a few top poses for rescoring , if the binding is non - specific or the binding mode is not well determined by the docking calculation .
the decomposition of the calculated binding free energy has demonstrated that the electrostatic energy changes upon binding in both gas phase ( eele ) and gb solvent ( ggb - ele ) are much less than those of the compounds in table 1 . on average
the van der waals energy changes ( evdw ) upon binding are more than those of the compounds in table 1 . for penfluridol , the evdw is the main factor to determine the binding affinity .
pose top6 has the largest evdw , which results in the largest binding affinity for all thirteen docking poses . as we pointed out in the previous section
, the electrostatic interaction plays a crucial role in the ligand binding of the antithrombin active site due to the large number of positively charge residues present in the active site .
compound l1c4 and ntp possess strong electrostatic interactions with antithrombin and can bind much tighter than penfluridol .
penfluridol has only two hydrogen bond donors , nitrogen and oxygen ( n1 and o1 ) .
hydrogen bonding analysis of penfluridol and antithrombin complex reveals that it has fewer hydrogen bond pairs and are less persistent during the md simulations . for the pose top1
, there are only two hydrogen bonds ( o1-hd22-nd2 and n1-hd22-nd2 ) formed during the md simulation between penfluridol and two asparagine residues ( asn61 , asn414 ) .
the hydrogen bond distances for o1-hd22-nd2 and n1-hd22-nd2 are 2.86 and 2.92  , respectively .
the persistence of hydrogen bond o1-hd22-nd2 is 16.62% and that of hydrogen bond n1-hd22-nd2 is 36.11% .
only one hydrogen bond exists in the initial docking pose , where hydrogen bond o1-hd22-nd2 is formed between penfluridol and asn61 .
after ~18 ns of simulation , hydrogen bond o1-hd22-nd2 is split while hydrogen bond n1-hd22-nd2 starts to form .
the docking score is relatively high in this pose because of a hydrophobic interaction between a phenyl ring and val375 and a - interaction between a phenyl ring and arg221 .
the - interaction is lost and two hydrogen bonds form between penfluridol and arg221 after 36 ns .
however , the persistence is 55.15% , which is the largest among the thirteen docking poses .
the hydrogen bond is absent from the initial docking pose but appears quick after ~4ns .
pose top4 has four hydrogen bonds form between penfluridol and lys283/ser277 during the course of the md simulation .
it is very interesting that the four hydrogen bonds appear and disappear one by one sequentially .
any time point during the md simulation , there is only one hydrogen bond present between penfluridol and antithrombin .
pose top6 has the largest rescoring binding affinity amongst the thirteen docking poses , which is largely due to its binding in the hydrophobic groove of antithrombin .
pose top7 , top8 , top9 , top10 , near1 , near2 , and near3 have 0 , 1 , 3 , 3 , 4 , 5 , and 5 , respectively .
the poses nearest to the active site ( near1 , near2 , and near3 ) tend to have more hydrogen bonds .
the number of hydrogen bonds for the penfluridol docking poses is no more than 5 .
the good binders ( e.g. l1c4 and ntp ) each have more than 25 hydrogen bonds .
navarro - fernandez and colleagues have also predicted the interacting residues by using the flexscreen scoring function . for compound l1c4 ,
they identified lys11 , asn45 , arg46 , arg47 , lys114 , and lys125 as key residues , but arg13 is missing from their list . for the compound l1c4 docking calculation , they used the receptor structure from the x - ray crystal structure of antithrombin complexed with compound ntp .
a few docking programs are able to have set side - chains of key residues in the receptor as flexible .
however , most docking programs can not sample the larger configuration space for the whole ligand - receptor complex . from our md simulations , we find that arg13 is quite flexible and can be adjusted to accommodate both large ( e.g. compound ntp ) and small ( e.g. compound l1c4 ) compounds , which shows the advantages of using md simulations over simple docking calculations .
judging from the hydrogen bond analysis on compound l1c4 and ntp , arg47 , arg13 , and asn45 play crucial roles in the antithrombin binding process .
antithrombin provides multiple sulfate / phosphate binding sites , consisting of mostly positively charged residues ( arginine , lysine ) and neutral charged residues that can provide rich hydrogen bond donors / acceptors ( asparagine ) .
all four phosphate groups of compound l1c4 form hydrogen bonds with antithrombin while not all the sulfate groups of compound ntp can form hydrogen bonds with antithrombin .
as pointed out above , introducing a positively charged group in the ligand will result in a penalty for solvation free energy .
if the positively charged group also can not form favorable interactions ( e.g. hydrogen bonding ) , ligand efficiency will be reduced , explaining why compound l1c4 has higher ligand efficiency than compound ntp . comparing the results from single - structure and ensemble - average mm / gbsa rescoring , the latter yields more accurate results .
the ensemble - average mm / gbsa rescoring ranks the binding affinities of antithrombin ligands in the order that agrees with the experimental results .
the advantage of ensemble - average mm / gbsa rescoring is that the binding affinity is averaged from an ensemble of structures extracted from long - time md simulations .
long - time md simulations can explore more configuration space and find energetically favorable configurations , which could offset the bias of initial structures .
arg13 is observed to form hydrogen bonds with the phosphate group of compound l1c4 after 8 ns of md simulation .
have also identified the importance of md sampling for recognizing the correct binding poses , although they carried out relatively short - time ( 300ps ) md simulations for the top three poses of 13 protein  ligand complexes and extracted only 60 snapshots from the trajectory for the ensemble - average mm / gbsa rescoring .
they concluded even short md simulation could help mm / gbsa to identify the experimentally determined binding conformations in most cases .
the ensemble - average mm / gbsa rescoring method is relatively accurate compared to single - structure mm / gbsa .
however , the method to obtain an ensemble of structures , long - time md simulations , is computationally intensive .
since the drug - like virtual libraries often contain millions of compounds , high - throughput virtual screening could be very costly , if using a more accurate , but expensive , method at the very beginning of the screening . to bridge the gap , our virtual screening pipeline uses a down - select scheme , a hierarchical approach , to screen large virtual compound libraries . a standard procedure to run the pipeline is to down - select compounds after they pass each screening method as implemented in the pipeline .
top ranked poses of down - selected compounds after docking are rescored using a single - structure mm / gbsa rescoring method . from this study , we know it is important to keep a few top docking poses for rescoring because the best binding pose in rescoring is often not the top docking pose , especially for non - specific binding or not well - determined binding modes in the docking calculation .
finally , the most expensive ensemble - average mm/ gbsa rescoring method in the pipeline can be applied to an amenable number of compounds down - selected after single - structure mm / gbsa rescoring , providing the accuracy needed for a fewer number of compounds .
in this article , we introduce a new addition , ensemble - average mm / gbsa rescoring , to our virtual screening pipeline . as a proof of concept , we calculated the binding affinities of eight antithrombin ligands by employing the previous single - structure mm / gbsa rescoring method and newly developed ensemble - average mm / gbsa rescoring method .
a total of 20 md simulations were performed for the ensemble - average mm / gbsa rescoring method .
the correlation coefficient of calculated and experimental binding affinities is improved from 0.36 to 0.69 when using ensemble - average mm / gbsa rescoring .
the rank order of calculated binding free energies using ensemble - average mm / gbsa rescoring exactly matches the experimentally derived free energies .
we demonstrate that long - time md trajectory can explore more configuration space and find energetically favorable configurations so that it can offset the bias of initial structures and improve the accuracy of binding affinity prediction .
the electrostatic interactions in both solute and solvent contribute favorably to the binding free energy . adding more negatively charged groups to the ligand provides more favorable electrostatic energy change .
the penalty can be compensated for by forming more hydrogen bonds as more negatively charged groups are added into the ligand .
the negatively charged groups added to the ligand must actively interact with the receptor by forming hydrogen bonds to achieve high ligand efficiency .
compound l1c4 has higher ligand efficiency because it uses all its phosphate groups to form hydrogen bonds with antithrombin while compound ntp does not .
a hydrophobic compound , penfluridol , has been docked to the antithrombin , which has abundant charge residues on the surface and in the active site of the protein .
the electrostatic energy changes , upon binding penfluridol in both gas phase and gb solvent , are much less than those of good binders .
the van der waals energy change upon binding is the dominant factor for the penfluridol binding affinity .
the number of hydrogen bonds present in the penfluridol md simulations are fewer than that of the good binders . we have also demonstrated that including a few top docking poses for rescoring is important especially for non - specific binding or binding modes that are not well - determined in the docking calculation .
we have previously reported 7 long - time md simulations as part of conference proceedings , and here we expand upon this work to disclose a total of 20 long - time simulations .
 background : a high - throughput virtual screening pipeline has been extended from single energetically minimized structure molecular mechanics / generalized born surface area ( mm / gbsa ) rescoring to ensemble - average mm / gbsa rescoring . 
 the correlation coefficient ( r2 ) of calculated and experimental binding free energies for a series of antithrombin ligands has been improved from 0.36 to 0.69 when switching from the single - structure mm / gbsa rescoring to ensemble - average one . 
 the electrostatic interactions in both solute and solvent are identified to play an important role in determining the binding free energy after the decomposition of the calculated binding free energy . 
 the increasing negative charge of the compounds provides a more favorable electrostatic energy change but creates a higher penalty for the solvation free energy . 
 such a penalty is compensated by the electrostatic energy change , which results in a better binding affinity . 
 a highly hydrophobic ligand is determined by the docking program to be a non - specific binder.results : our results have demonstrated that it is very important to keep a few top poses for rescoring , if the binding is non - specific or the binding mode is not well determined by the docking calculation . 
Background:
Results:
INTRODUCTION
METHODs
ReSULTS AND DISCUSSION
CONCLUSION
DISCLOSURE
CONFLICT OF INTEREST