molecular docking algorithms attempt to determine the binding modes of small organic molecules relative to a biomolecular receptor and to evaluate a score representing their relative binding propensity . in an effort to find novel binders for  hit  identification in structure - based drug discovery , virtual
screening ( vs ) entails utilizing a docking algorithm to rank large libraries of compounds .
receptor coordinates are most commonly provided by x - ray crystallography experiments as well as homology modeling or computer simulation .
the nature of the receptor model employed affects the predictive performance of dock - based approaches , as different conformations can produce alternative rankings of possibly active and inactive compounds , and only approximates the dynamic process occurring . although a number of strategies to incorporating protein flexibility have been developed in this context ( see e.g. refs ( 14 ) and references therein ) , defining protocols to select receptor structures for blind vs predictions is difficult . modeling the natural dynamics of a protein for ligand - binding events can benefit from methods that use multiple target configurations , so - called ensemble approaches , but not without limitations and trade - offs between sufficient model reliability and computational costs .
previous studies focused on crystallographic and homology models to examine single receptor effects on vs ranking performance as well as on the advantages of using multiple structures . yet ,
to our knowledge , no critical assessment of vs predictive power using individual receptor conformations from molecular dynamics ( md ) simulations has been reported to date .
this raises two general questions : are snapshots from md simulations predictive , and how do they compare to x - ray structures in affecting vs predictive power ? how do structures from the different types of md ensembles affect vs predictions ?
vs of md snapshots have been successfully used for pose prediction and compound library ranking . in some cases ,
clustering algorithms can alleviate computational costs by reducing the md ensemble with no significant loss of information for vs approaches .
for example , rare protein configurations have been shown to determine ligand binding in fkbp.(21 ) in other cases the dominant , frequent protein configurations are those promoting the best binding conditions for a variety of ligands . in the present study vs predictive power , using md snapshots and x - ray structures for two model systems , was explored .
the first model system selected was hiv-1 reverse transcriptase ( rt ; figure 1a ) .
rt catalyzes the transcription of the single - stranded rna viral genome into a double - stranded dna form and is essential for hiv replication . as a major drug target ,
rt is the subject of substantial structural biology efforts , resulting in more than a hundred related crystal structures to date .
together with computational studies , the heterogeneous properties of rt structures suggest considerable plasticity , which has been interpreted in the context of its function as both a dna polymerase and ribonuclease .
current fda - approved anti - rt drugs bind to one of two identified sites : the polymerase active site or a nearby hydrophobic allosteric site targeted by non - nucleoside reverse transcriptase inhibitors ( nnrtis).(24 ) the nnrti binding pocket ( nnibp ; figure 1c ) was the focus of the current work reported here , as it is of significant pharmaceutical interest and was suggested to be remarkably flexible , fluctuating between a  collapsed  inhibitor - free state and an  open 
inhibitor - bound state ( see , e.g. , refs ( 24 and 25 ) and references therein ) .
moreover , the nnibp has been shown to bind to a broad range of nnrtis , which bear structurally diverse scaffolds and was considered representative of allosteric binding sites.(24 ) protein receptors considered in this study : ( a ) rt and ( b ) w191 g overall representations on the same scale .
secondary structure elements and the location of the binding sites are highlighted ( red : helices ; cyan : sheets ; and gray : loops and turns ) .
insight views for : ( c ) the rt nnrti binding pocket ( nnibp ) with nevirapine bound and ( d ) the w191 g cation - binding pocket with 2a5mt bound .
ligands ( balls and sticks ) and pocket volumes ( blue surfaces ) are also shown .
the second model system was the less flexible w191 g artificial cavity mutant of cytochrome c peroxidase ( w191 g ; figure 1b ) for which ligand binding was characterized by x - ray crystallography and isothermal titration calorimetry experiments(28 ) as well as by computational techniques . taken as a whole , numerous x - ray structures and md simulations showed that no major backbone conformational rearrangement occurs upon binding of a large number of small cation compounds .
an exception to this general observation is made for three bulky ligands , which upon binding were found to induce the partial opening of the 190195 gating loop . in all cases ,
the electrostatic interaction between the cation ligands with asp235 , at the bottom of the binding pocket ( figure 1d ) , is crucial for binding.(31 ) in line with these distinctive properties , w191 g ligand binding is a prototype of lowest - possible change of receptor entropy.(31 ) rt and w191 g binding sites are ideal for evaluating the dependence of vs predictive power on diverse input receptor structures .
the subsequent sections describe analysis of blind vs to different ( unbound and bound ) md and crystallographic structures of the rt and w191 g systems .
predictive power was characterized using receiver operating characteristic ( roc ) curves to determine the probability of ranking actives over inactives . in this regard ,
the goal of this study was to determine if individual md snapshots were preferable for ranking actives over inactive compared to individual crystallographic structures .
the current article presents a procedure for nonarbitrary assessment and comparison of vs predictive power , of general interest for computational drug discovery .
these crystallographic conformations represent a combination of 10 diverse nnrti - bound and five nnrti - free states . in agreement with previous studies ,
the nnibp was defined as including pro 95 , leo 100 , lys 101 , lys 103 , val 106 , val 108 , val 179 , tyr 181 , tyr 188 , gly 190 , phe 227 , trp 229 , leu 234 , his 235 , pro 236 , and tyr 318 and is displayed in figure 1c . for w191 g , a total of 49 x - ray crystallographic structures were deposited in the protein data bank ( pdb ) and considered for our study ( 45 of which were cocrystallized with a ligand ; for a complete list of the corresponding pdb entries see supporting information , table si ) . out of the 45 binders reported
, modest binding site rearrangements are limited to amino acid side chains in 42 cases , thus excluding the peptide backbone . in agreement with previous studies ,
w191 g binding pocket includes his 175 , leu 177 , lys 179 , thr 180 , pro 190 , asn 195 , phe 202 , met 230 , met 231 , leu 232 , and asp 235 displayed in figure 1d .
pro 190 and asn 195 residues are the hinge points of w191 g cavity gating loop .
md simulations of four hiv-1 rt systems ( 983 residues ) in explicit water at 300 k were performed using the gromacs simulation software , the gromos 53a6 force field , and a compatible simple point charge ( spc ) water model.(39 ) the four different systems included two bound systems , 1vrt(40 ) structure with -apa bound to nnibp ( -apa bound ) and 1vrt structure with uc-781 bound to nnibp ( uc-781 bound ) as well as two unbound systems , 1vrt structure with the nnrti extracted from nnibp ( unbound open ) and 1dlo structure(41 ) ( unbound closed ) . for each system , four independent 30 ns trajectories were generated , with the final 12.5 ns of each used in our analysis .
md snapshots were extracted every 20 ps from this cumulative total of 50 ns of simulation time , yielding 2500 snapshots for each system and a total of 10  000 rt snapshots . for computational details
three md trajectories of the w191 g cavity mutant of cytochrome c peroxidase from escherichia coli ( w191 g ; 290 residues ) in explicit water at 300 k were generated using the gromos05 software,(42 ) the gromos 45a4 force field,(43 ) and a compatible spc water model.(44 ) the three systems included an initial closed - gate w191g-2a5mt complex from pdb : 1aen ( 2a5mt bound ) , an initial open - gate
w191 g from pdb : 1ryc after extracting ligand ( unbound open ) , and an initial closed - gate w191 g from pdb : 1aa4 ( unbound closed).(27 )
md snapshots were extracted every 20 ps , for a total of 2500 snapshots from 50 ns of cumulative trajectory for each complex , totaling 7500 w191 g receptor structures .
. for computational details refer to refs ( 23 and 30 ) . in the case of rt ,
20 diverse known inhibitors were combined with the national cancer institute diversity set ii ( http://dtp.cancer.gov ) , assumed to be decoys for the rt ligand set .
the w191 g ligand set comprised 53 known actives and 21 known inactives . in all cases , bond orders , stereochemistry , hydrogen atoms , and protonation states
were assigned by ligprep.(45 ) properties of the ligand sets are summarized in table 1 .
these compound sets were chosen to represent prototypes of challenging sets for blind vs. see supporting information for a complete list of x - ray structures , actives , and decoys considered .
for all receptor structures , crystallographic and trajectory water molecules , ions and ligand compounds were removed .
this blind setup reflects the typical scenario for vs of newly discovered protein targets , in which no additional chemical information can be used a priori for docking .
receptors were processed using the default protein preparation workflow , which employs a restrained , partial energy minimization .
the glide sp algorithm was employed using a grid box volume of 10  10  10  .
all structures were fitted using the c atoms of the binding pocket , as defined above , and grid center coordinates were based on this alignment . only the best , or lowest energy , pose for each docking run was retained . the use of glide as a docking software might influence the quantitative results of this study , however , the more general qualitative conclusions are likely transferable to different algorithms and scoring functions .
a total of 2500 md snapshots for each of the 7 systems were processed as well as 15 rt and 49 w191 g x - ray structures , totaling 17564 different receptor docking runs performed and analyzed in this study .
after the algorithm ranks the compounds , the top x will be pursued further , by experimental assay or as a starting point for rescoring calculations , the rest are discarded .
the top x are thus deemed positives , and the rest negatives . if the activity of the compounds ranked is previously known , then such classifications are either true or false . using the four categories true positive ( tp ) , false positive ( fp ) , true negative ( tn ) , and false negative ( fn ) , the true positive rate ( i.e. , the selectivity):and the false positive rate ( i.e. , 1-specificity):are determined based on the x threshold .
the roc plots eq 1 with respect to eq 2 as x is perturbed , producing a curve that is threshold independent .
the integral of a roc plot , or the area under the curve ( auc ) , gives the probability that a classifier , in this case the docking algorithm , will rank a randomly chosen active over a randomly chosen inactive .
figure 2 summarizes the practical aspects of comparing vs predictive power for two conformations of a receptor .
different a and b receptor conformations used in docking affect compound ranking , thus prediction selectivity and specificity .
the integral of a roc curves ( also named auc ) is the metric employed for nonarbitrary comparison of vs predictive power .
the mean auc for different groups of structures , bound and unbound , was compared for quantitative assessment .
the best , worst , and mean auc for 10 and 45 bound x - ray structures as well as 5 and 4 unbound x - ray structures for rt and w191 g , respectively , were calculated for comparison . all of these x - ray groups were compared to averages based on 2500 snapshots for each of the 7 systems .
docking the compound sets to 2500 md trajectory snapshots results in a distribution of auc values for each of the seven systems .
these distributions arise from the diversity of the conformational space sampled by the receptor and allow for quantification of bound versus unbound vs predictive power . in addition
, different protein systems auc distributions can be compared on a qualitative basis , keeping in mind that both receptor configuration sampling and ligand diversities determine such representations .
auc values of 0.5 , corresponding to a random ranking of the compounds , or smaller signify worse than random predictive power for a give ligand set .
figure 3 shows auc distributions based on four rt and three w191 g trajectories and compared with the range of and average corresponding auc values from individual x - ray bound and unbound structures .
all trajectory distributions are largely unimodal with the exception of the w191 g unbound , started from the pdb structure 1ryc .
the presence of a second peak can be explained considering the comparatively greater flexibility of the open 190195 gating loop.(30 ) normalized distributions of auc values for blind vs of rt ( left panels ) and w191 g ( right panels ) systems based on seven molecular dynamics ensembles .
auc value ranges ( delimited by vertical black lines corresponding to worst and best values ) and the auc averages ( vertical red lines , n is the number of crystallographic structures available ) from corresponding x - ray structures are also displayed as reference .
area under the distribution to the right of the vertical red lines indicates md structures with better predictive power than average x - ray structures from corresponding ensembles .
a pronounced shift toward lower auc values and poorer predictive power can be observed by comparing the bound systems peaks of 0.76 and 0.78 auc for -apa and uc-781 bound systems , respectively , and the unbound ensemble peak of 0.43 and 0.44 auc for open and closed unbound simulations , respectively .
this is consistent with previous studies suggesting improved vs predictive power of complex structures compared with unbound receptors.(15 ) in the case of w191 g , distributions for bound and unbound systems are more similar , consistent with the fact that receptor relaxation upon ligand binding is less pronounced for more rigid binding sites .
the mean , maximum , and minimum auc values from the bound and unbound ensembles are reported in table 2 .
interestingly , in all cases , the best md auc value outperforms the best x - ray available .
however , a significant percentage of the md aucs was less predictive than the worst x - ray structure available .
the mean values indicate that on average the predictive power of md snapshots is poorer than that of the x - ray structures .
an exception to this observation is made for rt unbound , for which the mean values are comparable .
refer to the methods section for details on the different ensembles considered and to figure 3 for corresponding distributions .
the range of predictive power of x - ray structures can also be seen in figure 3 with black vertical lines at maximum and minimum auc values found in this study .
the bound x - ray ranges for both systems are wider than unbound x - ray auc ranges , although rt has a considerable larger range , consistent with greater flexibility in this case .
interestingly , similar distribution ranges of ca . 0.5 auc values for md snapshots are found for both systems .
the predictive power of md versus x - ray structures is represented in figure 4 in terms of the relative number of md snapshots with auc values larger than the average value from corresponding x - ray models .
this is a measure of the advantage of using md structures versus using a randomly selected static x - ray model . for rt ligand - bound systems at least 23% of md snapshots
lead to improved vs prediction , this value being somewhat higher for the -apa versus uc-781 system ( 31% ) . for rt unbound systems , at least 50% of md snapshots outperform static x - ray models with similar values for unbound open ( 51% ) and unbound closed ( 50% ) simulations .
relative number of molecular dynamics snapshots with higher predictive power ( auc ) than the mean auc from all available corresponding x - ray structures .
see also supporting information , table sii , for detailed data . for the rigid w191 g ,
the majority of the md snapshots were less predictive than the average x - ray auc , with 12% for the bound simulation better than the average x - ray structure auc , and less than 1% for both unbound simulations compared to the average auc of bound x - ray structures .
it is important to stress that auc values are randomly distributed along the seven trajectories ( not shown ) .
therefore , there is no correlation of auc value with the evolution of the md run from the same initial x - ray coordinates .
as described previously , the auc value of 0.5 corresponds to the predictive power of randomly ranking compounds .
figure 5 summarizes the vs predictive power of md snapshots in terms of the number of md snapshots with auc values greater than 0.5 or random .
rt percentages of predictive conformations were 99.8% -apa and 99.9% for uc-781 ligands versus 25.3% and 20.9% for unbound open and closed , respectively .
w191 g percentages of predictive conformations were 94.0% for the bound system versus 71.6% and 88.2% for unbound open and closed , respectively .
for both systems , at least 94% of bound md snapshots were more predictive than random .
the unbound simulation predictive structures were less in number , particularly for the unbound rt structures .
relative number of molecular dynamics snapshots with higher predictive power ( auc ) than random ranking .
remarkably , in both unbound rt simulations , approximately 50% of the structures processed were more predictive than the average of the unbound x - ray aucs .
more significantly , 21% of the closed unbound and 25% of the open unbound simulations were more predictive than random , where as all of the x - ray unbound structures were worse than random .
this case indicates that md can move a conformation seemingly not amenable to docking , into the predictive range .
mcgovern et al . found that the majority of x - ray structures studied identified known binders over nonbinders better than random selection of rank.(15 ) the x - ray structures for the two systems studied here were all more predictive than random , having a auc of greater than 0.5 , with the exception of unbound rt structures which were all worse than random .
the reason for the poor performance of unbound rt x - rays is likely due to the occluded binding site , where key side chains must move out of the way in order for a ligand to fit .
figure 6 reports a comparative analysis of backbone fluctuations for rt and w191 g systems .
however , the all - atom binding site fluctuations are comparable for both systems ( supporting information , figure s1 ) , and a pronounced decrease in side chain flexibility upon ligand binding can be observed , when comparing bound and unbound trajectories . root - mean - square fluctuations ( rmsf ) of the backbone c atom positions from the seven molecular dynamics simulations .
rt ( left ) and w191 g ( right ) , while ( a and b ) show bound simulation fluctuations and ( c and d ) show unbound simulation fluctuations .
rmsf values were calculated using sampling of 50 ns for each of the seven cases . for analysis of these simulations
see ref 25 ( rt ) and refs 23 , 30 ( w191 g ) .
looking at fluctuations of all binding site atoms and the -carbon subset represents motion with respect to the entire binding site and the backbone chain that makes up the area .
figure 7 plots the average binding site fluctuations for each of the trajectories , comparing it with the average auc of the ensemble .
a correlation was found between the average predictive power < auc > and the averaged binding site flexibility < rmsd>. averaged fluctuations and predictive power values ( supporting information , table siii ) show generally that the more flexible the binding site , the less predictive the md ensemble .
this is in agreement with previous findings that bound or holo structures , traditionally more stable and less flexible , are more amenable to predictive docking experiments than unbound or apo structures.(15 ) relationship between ensemble averaged predictive power , < auc > , and ensemble averaged binding site flexibility , < rmsf>. values were calculated using sampling of 50 ns for each of the seven cases .
see supporting information , table siii , for data values . while the two systems display variation in flexibility , the different chemical nature of the compounds docked as well as of the active sites should also be noted .
compounds in the rt set represent a more diverse range of physicochemical properties than for w191 g .
additionally , w191 g ligand binding is dominated by a specific electrostatic interaction , with all known binders carrying a positive charge .
while ligand ranking might be equally challenging for both systems , the more pronounced difference between ligand and decoys for w191 g may in part explain a higher overall vs predictive power in this case .
analyses of root - mean - square deviation ( rmsd ) from starting structure , binding site volume , radius of gyration , cognate ligand size , mean score of top 1% of docked compounds , rmsd binding site clustering , and flexibility descriptors from principal component analysis were not able to distinguish better performing snapshots from weakly performing snapshots .
overall , it was not possible to identify any evident relationship between characterization of md snapshots and vs predictive performance .
the results suggest that md conformations can improve vs results compared to the exclusive use of x - ray structures , yet a general system independent protocol to this goal is still challenging .
the present data strongly confirms previous observations that the vs  quality  of a protein receptor configuration is strictly dependent on the chemical properties of the system and the ligand sets considered .
the only general consideration demonstrated by this study is that no simple property can be used to separate the best receptor structures from the ensemble , and as such , none should be discarded a priori . while no method for selecting md structures before hand was determined , exploratory screening using the auc measurement on an md ensemble will be used to identify the best md conformations prior to a more extensive vs computation , in a hierarchical fashion in future work .
additionally , while individual receptor structures were assessed here , evaluation of efficient methods for combining information from multiple ensembles is on going .
one practical concern for users seeking enrichment is that generating more conformations using md simulations may be computationally time intensive .
as technology advances , time scales achieved in even one days worth of simulation time are rapidly increasing , and it is plausible that in the future generating and choosing a receptor structure may be as practical as choosing and generating the number of compounds in a vs database ; this was the main motivation for evaluating the different ensembles , and indeed it was shown here that md can provide snapshots that are different from crystal structures . in most cases , the extent of md sampling will determine the diversity of conformations , relying on how well the sampling has explored phase space .
enhanced sampling methods , such as accelerated molecular dynamics , allow for more rapid exploration of conformational space by reducing energy landscape barriers and subsequently effective simulation time scale.(48 ) finally , one could argue that the user becomes more intimately familiar with the target through running a dynamics simulation ; standard practice for docking protocol often involves manual review and processing as well as familiarity with the system .
this was highlighted by head s results from sampl-1 cross - docking challenge , which showed that successful use of docking programs requires human expertise and experience.(49 ) furthermore , some publications of md vs studies have appeared subsequent to an md study that aimed to biophysically characterize a target ; in such cases simulation data were available for vs , and the computation time is dependent just on the analysis or selection of structures , which is done in the case of multiple x - ray structures as well .
this study presented an explorative analysis of blind virtual screening predictive power for ligand bound and unbound molecular dynamics snapshots of two systems with diverse molecular flexibility and binding site properties .
predictive power was quantified based on virtual screening using receiver operating characteristic curves to determine the probability of ranking actives over inactives .
results from predictive power analysis were compared with crystallographic structure results as well as with respect to random compound rank .
most interestingly , the best molecular dynamics snapshot was systematically more predictive than the best x - ray structure in all cases examined .
however , particularly for the less flexible of the two systems considered , a significant number of molecular dynamics snapshots led to less predictive virtual screening results than the least predictive x - ray structure available .
therefore , the advantage of molecular dynamics simulations for virtual screening seems to rely on the enrichment that a reduced number of receptor structures supply , rather than the whole configurational ensemble generated . while no method for selecting the best molecular dynamics structures a priori was found , exploratory screening will be used to identify the best conformations in a hierarchical approach . for the more flexible rt system ,
enhanced sampling methods are being employed to improve virtual screening predictive power to an even greater extent .
 molecular dynamics ( md ) simulation is a well - established method for understanding protein dynamics 
 . conformations from unrestrained md simulations have yet to be assessed for blind virtual screening ( vs ) by docking . 
 this study presents a critical analysis of the predictive power of md snapshots to this regard , evaluating two well - characterized systems of varying flexibility in ligand - bound and unbound configurations . 
 results from such vs predictions are discussed with respect to experimentally determined structures . in all cases , 
 md simulations provide snapshots that improve vs predictive power over known crystal structures , possibly due to sampling more relevant receptor conformations . additionally , md can move conformations previously not amenable to docking into the predictive range . 
Introduction
Materials and Methods
Results
Discussion
Conclusion