docking is a well - established computational 
 technique often used 
 at the early stages of structure - based drug discovery . in general , the purpose of a docking experiment is to predict the 
 optimal binding geometry ( pose ) of a ligand ( typically 
 a small organic molecule ) with respect to a receptor ( typically a protein drug target ) , and to provide an estimate of 
 the interaction affinity relative to other ligands . dock is one such program designed for binding pose prediction 
 and , when used in a virtual screening capacity to rank - order large 
 databases of candidate molecules , is a powerful approach for lead 
 discovery . numerous dock successes have been reported ; recent examples include identification of leads targeting ctx - m 
 -lactamase , 2-adrenergic receptor , xpa protein , and the riboswitch of mrna , among others . examples of successful virtual screens from 
 our laboratory employing the most recent version of dock ( version 
 6 ) include the drug targets hivgp41 and fatty acid binding protein . other related programs , including autodock , autodock vina , gold , glide , surflex , and 
 flexx , are widely used for similar goals 
 and have documented analogous successes ( reviewed in yuriev et al . ) . root - mean - square deviation ( rmsd ) , a 
 measure of intermolecular 
 differences in position and conformation , has broad applications in 
 computational biology , including testing the efficacy of docking programs . for example , in a typical docking validation experiment , a cocrystallized 
 ligand is removed from its binding site , then multiple candidate binding 
 poses are predicted using a specific protocol . the positive - control 
 experiment is successful if the heavy - atom rmsd of the top - ranked 
 pose is within a certain tolerance from the original ligand 
 position ; a typical cutoff for success is 2.0  . however , if 
 a molecule contains symmetric functional groups or is itself symmetric , 
 inflated rmsd values and a false negative result can occur . as an example , docking 1,2-dichlorobenzene to 
 t4 lysozyme ( pdb code 2oty ) can yield chemically 
 equivalent binding poses , yet different rmsd values ( figure 1 ) . in one scenario , the pose rmsd is less than 2.0 
  from the crystallographic reference and is considered a docking 
 success . in the second scenario , the pose is inverted over its axis 
 of symmetry , yielding an rmsd greater than 2.0  , which is reported 
 as a docking failure . dock6 was recently reported to successfully 
 redock between 63.8% ( over 780 systems ) and 66.3% ( over 147 systems ) of known 
 crystallographic ligand receptor complexes in a series of flexible 
 ligand docking experiments using a standard energy - based scoring function . 
 a subset of the failures from these experiments could be attributed 
 to molecular symmetry in the ligand . when symmetry is accounted for , 
 an increase in success of approximately 48% these trends demonstrate 
 that neglecting to account for symmetry in ligands can negatively 
 impact the accurate evaluation of new docking programs , protocols , 
 search methods , or scoring functions . docked poses , shown in orange , are inverted 
 approximately 180 about the axis of symmetry . , there remains a 
 surprising lack of published details regarding specific functions 
 used for symmetry correction , as was recently noted by brozell et 
 al . a notable exception is the publication 
 describing autodock vina , which reports 
 the exact procedure used to compute what the authors refer to as the minimum - distance rmsd . however , while functionally simple 
 to implement , minimum - distance rmsd in some cases may not enforce 
 a one - to - once atom correspondence . under these circumstances , some 
 atoms could be used multiple times for the rmsd calculation while 
 others could be neglected , which may not be desirable ( see theoretical methods ) . as an alternative approach , 
 our group has implemented a symmetry - correction strategy into the 
 dock6 program based on the method described 
 by kuhn and later by munkres known as the hungarian algorithm . a similar strategy was recently reported by 
 helmich et al . to find the best overlap 
 in alignment between two conformations of the same molecule , and by 
 ellingson et al . to efficiently maximize 
 the number of atom superpositions in an overlay of protein binding 
 sites . in this report , we demonstrate that our implementation 
 of the hungarian 
 algorithm effectively corrects for symmetry in molecular docking which 
 provides the means to assess the impact of various algorithmic or 
 protocol changes as the dock program continues to evolve . further , 
 we discuss an extension of the algorithm for adding chemical diversity 
 to a pool of molecules that are assembled in a de novo design version 
 of dock currently under development , enabling improved sampling of 
 diverse chemotypes . additional applications presented include rank - ordering 
 compounds relative to a known reference to aid virtual screening and 
 computing symmetry - corrected rmsds to analyze molecular dynamics simulations . the standard pairwise rmsd function , 
 as it is most commonly used in structural biology , quantifies the 
 difference between two poses of the same molecule . in this context , 
 differences between poses can include rigid geometric rotations and 
 translations within the six standard degrees of freedom ( three rotational 
 and three translational ) , plus any changes in internal degrees of 
 freedom ( such as dihedral rotations ) . the standard rmsd function ( rmsdstd ) is defined as follows:1where molecule poses a and b are sets comprising individual 
 atoms : a = { a1,a2 , ... , an } and b = { b1 , b2 , ... , bn}. atoms ai and bi can be 
 interpreted as cartesian coordinates 
 or positional vectors : for example , ai = { xa , i , ya , i , za , i}. the total number of heavy 
 atoms in the molecule is n. historically , this 
 function is dependent on a predefined one - to - one correspondence between 
 atoms in two poses ( a and b ) of 
 the same molecule . for practical purposes , this is typically achieved 
 by matching atom indices or evaluating the sets of atoms in the same 
 order . in the case of highly symmetric molecules , however , it would 
 be advantageous to impose a dynamic atom - matching algorithm prior 
 to the standard rmsd calculation that is not subject to the constraints 
 of arbitrary atom numbers or orders , such that physically indistinguishable 
 atoms would be matched . in an attempt 
 to overcome this 
 limitation , trott and olson recently 
 described a modified function termed minimum - distance rmsd ( rmsdmin ) , computed as follows:2where3 in this method , 
 atoms ai from 
 reference pose a are iteratively compared to all 
 atoms bj from pose b that are of the same element 
 type . the minimum distance identified defines the atom correspondence , 
 which is subsequently used in the root - mean - square calculation . the 
 same procedure is repeated using pose b as the reference , 
 and the maximum of the two calculations is reported as the rmsdmin . however , as this function does not enforce 
 a one - to - one correspondence , some atoms may be used more than once 
 in the calculation while others that are far removed in coordinate 
 space have the potential to be neglected . therefore , while this method 
 is effective , it can under - represent the true rmsd . in order 
 to account for all atoms in each pose a and b , the minimum - distance rmsd could be modified 
 to enforce a one - to - one correspondence through removal of atom pairs 
 after they have been matched via following the steps : ( 1 ) for each i , compute the minimum distance : minj||ai  bj|| , ( 2 ) remove ai and bj from 
 atom sets a and b , respectively , 
 and ( 3 ) repeat steps 12 until no atoms remain to be matched . once the complete atom correspondence is determined , the rmsd would 
 be computed similar to eq 3 . however , while 
 this modified approach would ensure that no atoms are neglected from 
 the final rmsd calculation , it is dependent on the order in which 
 atoms are evaluated and may give a nonoptimal solution . as an alternative approach , 
 the optimal one - to - one correspondence between sets 
 of atoms can be solved using the hungarian algorithm , a method for solving the minimum assignment problem . in this approach , every atom from pose a will be matched to exactly one atom from pose b in such a way that minimizes the sum of distances between 
 all atom pairs , and no outlying atoms will be neglected from the rmsd 
 calculation . the optimal - correspondence rmsd ( rmsdcor ) , 
 termed here the symmetry - corrected rmsd , can be defined as follows:4where atom ai is matched uniquely to 
 an atom bj that was predetermined by solving the optimal correspondence 
 corj||ai  bj|| . briefly described , 
 the hungarian algorithm solves this correspondence in four steps : 
 ( 1 ) form a cost matrix m where each matrix element mij is the squared norm of atom 
 positional vectors ai and bj . subtract the lowest value in each 
 row from all elements in the same row and the lowest value in each 
 column from all elements in the same column . ( 2 ) find the minimum 
 number of horizontal and vertical lines that can cover all zeroes 
 in the matrix . if the number of lines is equal to the size of the 
 matrix in one dimension , a set of independent zero - values exists , 
 and the positions of these zeroes represent the optimal - correspondence 
 solution ( 3 ) determine the smallest 
 matrix element that is not covered with a line . subtract that value 
 from all uncovered matrix elements and add that same value to all 
 matrix elements that are covered by two lines . ( 4 ) repeat steps 23 
 until the matrix is solved as described in step 2 . each iteration 
 of steps 23 constitutes a cycle of the algorithm . 
 a more thorough description and proof are documented by munkres . a very fast , o(n ) version of the 
 hungarian algorithm was implemented into 
 the source code of dock6 . disjoint distance matrices ( termed submatrices ) 
 are assembled for each atom type ( using tripos sybyl conventions ) in a given molecule , thus only atoms of the 
 same type can correspond to one another in the rmsd calculation ( see 
 figure 2 ) . all tests were performed on a dell 
 poweredge c6100 linux cluster consisting of intel x5660 2.8 ghz hex - core 
 nehalem - based processors . in this example , a theoretical molecule that 
 contains only four 
 sybyl atom types ( o.3 , f , c.3 , and n.3 ) is described . the search 
 space was reduced from all possible atom comparisons ( left ) to comparisons 
 between only atoms of the same atom type ( right ) . the amount of calculation 
 required , represented by blue shading , is significantly decreased 
 by dividing the molecule into subsets , or submatrices . pose - reproduction 
 docking 
 experiments were performed over a large test set using the flexible 
 ligand protocol described by mukherjee et al . the test set used is an expanded version of the sb2010 test set 
 ( n = 780 ) recently increased 
 to include 1043 receptor  ligand complexes , available for download 
 as sb2012 at www.rizzolab.org . key aspects of the test 
 set preparation and flexible ligand docking protocol include the following : 
 ligands were protonated within the program moe , then semiempirical am1-bcc charges were added with antechamber , part of the ambertools 
 package . receptors were protonated and 
 assigned force field parameters with tleap ( ambertools ) , 
 then the receptor structures were subjected to a short energy minimization 
 using sander ( amber11 ) with strong restraints on the heavy atoms . following standard 
 dock6 practices , the receptor surface was generated with dms , spheres were generated using the utility sphgen , and a box was defined that surrounded the spheres 
 plus an 8.0  margin in all directions . a docking grid was generated within the boundaries of the box 
 at 0.3  resolution with 69 van der waals exponents and 
 a distance - dependent dielectric of  = 4r . 
 during docking , a maximum of 1000 anchor orientations ligand positions 
 were minimized during growth with a simplex minimizer for a maximum 
 of 500 iterations or until the change in energy score between steps 
 was less than 1.0 kcal / mol . a maximum of 5000 fully grown conformers 
 were kept for each ligand , following best - first clustering with a 
 2.0  standard rmsd threshold . the resulting ensembles of candidate 
 poses ( 88  853 total poses ; average 85.2 poses per system ) were 
 evaluated using both the standard ( eq 1 ) and 
 symmetry - corrected ( eq 4 ) rmsd algorithms described 
 above . all rmsd values reported were measured between the docked pose 
 and the crystallographic position of the same ligand . rmsd values 
 are computed using heavy - atoms only  hydrogen atom positions 
 were not considered . as an alternative 
 application of the hungarian algorithm , we have adapted the code to 
 enable clustering and pruning between molecules of differing chemical 
 composition in order to aid in de novo design . depending on several 
 factors , the number of candidate molecules during de novo assembly 
 can quickly expand exponentially , thus rigorous pruning methods are 
 required to keep the search space tractable . to evaluate the impact 
 of the new pruning heuristic , experiments to reconstruct known molecules 
 and to construct diverse ensembles of molecules were performed using 
 a de novo version of dock currently under development ( further discussed 
 in balius et al . ) . prior to the 
 de novo experiments , two different groups of molecule fragment libraries 
 were prepared . first , all ligands from the sb2012 test set with exactly 
 seven rotatable bonds ( n = 103 ) were independently 
 disassembled on those bonds into 103 small , restricted fragment libraries . second , a subset of the zinc database consisting of approximately 1  10 representative drug - like molecules was used to form a larger , generic fragment library . new molecules were then assembled 
 de novo in the binding sites of each receptor from the sb2012 subset 
 using either ( 1 ) the smaller , restricted fragment libraries associated 
 with each system or ( 2 ) the larger , generic library from zinc seeded 
 with fragments from the restricted libraries . following each layer 
 of assembly , molecule ensembles were rank - ordered by the value of 
 their dock grid energy score , then pruned 
 using a best - first clustering method and the adaptation of the hungarian 
 algorithm . the criteria for pruning were as follows : all possible 
 atom matches ( using sybyl conventions ) between the two molecules were 
 identified , the optimal - correspondence of those matches was determined , 
 and the symmetry - corrected rmsd for matched atoms was computed ( hereafter 
 referred to as rmsd - matched ) . in addition , the atoms 
 that remained unmatched in either molecule were enumerated ( hereafter 
 referred to as # -unmatched ) . two molecules were deemed 
 to be sufficiently similar in cartesian space and chemical space if the rmsd - matched term was less than a user - defined 
 cutoff ( typically 2.0  ) , and the # -unmatched term was less than 
 a second user - defined cutoff ( typically 010 atoms ) . under 
 these conditions , the worse - scoring molecule would be pruned . in practice , 
 as the threshold for number of unmatched atoms increases , molecules 
 would be pruned more frequently . for clarity tanimoto coefficients were computed using a fingerprinting method 
 implemented into dock6 , inspired by the molprint 2d algorithm described 
 by bender et al . comparison between a gray reference molecule 
 ( same in all four 
 panels ) and four example molecules constructed by de novo design with : 
 ( a ) slight differences in functional group positions , ( b ) larger differences 
 in functional group positions , ( c ) slight differences in functional 
 group composition , and ( d ) larger differences in functional group 
 composition . in total , 
 88  853 docked 
 poses derived from the sb2012 test set were compared against 1043 
 crystallographic reference conformations using the standard and symmetry - corrected 
 rmsd functions . the present implementation of the hungarian algorithm 
 was made significantly more efficient through a process where submatrices 
 are prepared containing only atoms of the same sybyl atom type for 
 matching ( see figure 2 ) . in terms of the chemistry 
 of matching atoms , it follows that only atoms of an identical topology 
 should be matched to correctly account for molecular symmetry . in 
 order to verify that the hungarian algorithm could successfully solve 
 each minimum assignment problem in the data set , we measured the number 
 of algorithm cycles required for each calculation . the 88  853 
 poses were divided into 675  056 submatrices based on atom type 
 ( average 7.60 submatrices , or atom types , per molecule ) and the minimum 
 assignment was determined . scatter plot of the number of cycles of the hungarian 
 algorithm 
 ( y - axis ) required to solve the minimum assignment 
 problem for matrices of a given size ( x - axis ) . each 
 red point represents one submatrix ( n = 675  056 ) . the maximum number of algorithm cycles required to solve 
 the most 
 difficult matrix was 265 cycles , which occurred in a matrix of 33 
  33 elements . in addition , the largest matrix evaluated in the 
 data set was 35  35 elements , derived from a molecule with 138 
 total atoms , of which 35 atoms were the same sybyl atom type ( sp3 
 carbon , pdb code 1kqy ) . on average , the complexity of solving the minimum assignment 
 problem , as measured by the number of algorithm cycles required , increases 
 with increasing number of matrix elements . conceptually this is consistent 
 with increasing opportunity for alternative atom correspondence among 
 larger ligands . figure 4 also demonstrates 
 that all distance matrices from a very large data set were solvable , 
 exhibiting the robustness of the algorithm . finally , it is noteworthy 
 that rescoring all 88  853 poses ( 675  056 total submatrices ) 
 required 77  981 ms total of computer time . this value factors 
 to an average of 0.88 ms per pose or 0.12 ms per submatrix to solve 
 the minimum assignment problem . submillisecond time scales for evaluating 
 the symmetry - corrected rmsd are negligibly small when compared to 
 an average on - the - fly flexible ligand docking experiment in the program 
 dock6 , which under the present conditions required approximately 4.9 
 min per ligand . this millisecond time scale becomes an important consideration 
 when , as discussed below , the algorithm is adapted to enable pairwise 
 comparisons for de novo design . by design 
 of the algorithm , the symmetry - corrected rmsd can only be less than 
 or equal to the standard rmsd . that is , the atom correspondence will 
 not change from the standard correspondence if any change increases 
 the sum of distances between all atom pairs . a heat map of rmsd comparisons 
 for all 88  853 poses is shown in figure 5 . data here are grouped in 0.2  bins , red indicates areas of 
 high population , green indicates areas of low population , and white 
 areas are unpopulated . the heat map can be divided into four quadrants : 
 poses in quadrant i have both a low standard and symmetry - corrected 
 rmsd ( < 2.0  ) , thus would be considered successfully docked 
 poses by either metric . among the results , poses in quadrant ii have 
 high standard rmsd values ( 2.0  ) , which would ordinarily 
 be considered a docking failure , but they have low symmetry - corrected 
 rmsd values ( < 2.0  ) . thus , they are false - negatives , or poses that are rescued from failure by the 
 symmetry - correction algorithm . poses in quadrant iii have a high rmsd 
 by both metrics ( 2.0  ) , and are considered to be the true negatives . finally , quadrant iv ( false positives ) can not 
 contain any poses and should always remain empty . populations are binned in 0.2  increments , and are colored 
 from green ( low population ) to red ( high population ) as shown on the 
 right - hand side bar . dashed lines at x = 2  
 and y = 2  delineate four quadrants : i = true 
 positives , ii = false negatives , iii = true negatives , iv = false 
 positives . from figure 5 , it is evident that no poses 
 are assigned a symmetry - corrected rmsd greater than the standard rmsd , 
 consistent with expectations . the vast majority of poses fall near 
 to the diagonal , as indicated by the thick red band across the heat 
 map . however , while the average difference in rmsd between the two 
 metrics ( rmsdstd  rmsdcor ) is only 0.96 
  , 11  914 ( 13.4% ) poses differ by more than 2.0  , 
 and 5490 ( 6.2% ) poses differ by more than 3.0  . the poses that 
 differ by greater than 3.0  represent the most interesting cases 
 of highly symmetric molecules that the hungarian algorithm was designed 
 to uncover , as discussed in the next section . the vast majority of 
 poses in figure 5 , encompassed within the breadth 
 of the red band , are likely indicative of commonly occurring , smaller 
 corrections to symmetric functional groups . for example , as a point 
 of comparison , in a hypothetical molecule of 25 heavy atoms , a 180 
 rotation of a single phenyl group would result in relatively small 
 rmsd of 0.96  , while the 180 rotation of two phenyl groups 
 would result in a larger rmsd of 1.36  . to facilitate discussions 
 of differences between standard and symmetry - corrected rmsd , it is 
 convenient to analyze docking outcomes grouped by the four quadrants 
 as shown in table 1 . the results here are further 
 subdivided to specifically look at outcomes wherein only the top - scoring 
 pose was retained for each system ( n = 1043 ) , or 
 all poses that were generated were retained ( n = 
 88  853 ) . it is important to note that the poses analyzed here 
 follow a clustering protocol such that only a diverse set of poses 
 from each ligand docked to its target will be retained . these poses 
 are diverse in terms of both conformation and orientation within the 
 binding pocket . therefore , given the constraints of diversity , typically 
 only one pose would be expected to be the correct pose , or within 
 2.0  from the crystallographic reference . in addition , it is 
 expected that a small subset of poses , those that overlay within 2.0  
 from the reference but inverted over an axis of symmetry , will be 
 rescued by the symmetry - correction algorithm . all remaining poses 
 should be perceived as incorrect , or not within 2.0  from the 
 reference . from the perspective of the dock6 scoring function , the 
 correct poses and the rescued poses are likely to have indistinguishable 
 energy scores , and it is a statistical coin - flip that ultimately determines 
 which one is the top - ranked pose . among the data presented in 
 table 1 , these 
 expectations appear to hold true . specifically , approximately one 
 pose per system , or 1272 ( 1.4% ) , is identified as a true positive 
 ( quadrant i ) . the slight discrepancy between true positive poses ( n = 1272 ) and number of systems ( n = 1043 ) 
 can be attributed to the dock6 sampling algorithm , which may generate 
 two poses that are each less than 2.0  from the crystallographic 
 reference ligand , but greater than 2.0  from each other , in 
 which case both poses would be kept . in addition , 853 poses ( 1.0% ) 
 are identified as false negatives and rescued by the symmetry - correction 
 algorithm ( quadrant ii ) , which is consistent with the hypothesis that 
 only a small subset of poses can overlay with the crystallographic 
 reference , but be inverted over some axis of symmetry . finally , the 
 vast majority , or 86  728 ( 97.6% ) , of all poses are true negatives 
 ( quadrant iii ) . this demonstrates that the ensembles of candidate 
 poses generated by dock6 fully explore the boundaries of the binding 
 pockets and vary largely in conformation and orientation from the 
 crystallographic ligands . as expected , no poses are considered false 
 positives given that symmetry - corrected rmsd can never yield a higher 
 value than standard rmsd ( quadrant iv ) . the need to account 
 for symmetry becomes more important when only 
 the top - scored pose for each system ( n = 1043 ) is 
 examined , as is typical for most applications of docking and virtual 
 screening . while the top pose in table 1 is 
 ranked successfully by both rmsd metrics 67.7% of the time , which 
 is consistent with our laboratory s earlier statistics , the number of poses rescued through symmetry correction increases 
 from 1.0% to 5.0% . therefore , the true pose - reproduction success rate 
 using the flexible ligand docking protocols presented here becomes 
 72.7% ( 67.7% + 5.0% ) across all 1043 systems in the sb2012 test set . a scatter plot comparing standard and symmetry - corrected rmsd values 
 for the top pose from each system is shown in figure 6 . again , most of the top - scored poses fall near to the diagonal . focusing on quadrants i and ii ( figure 6 , insert ) , 
 specific examples of top - scored poses begin to emerge as individual 
 red points . the 20 most interesting cases ( out of 52 in quadrant ii ) 
 with large standard rmsd and small symmetry - corrected rmsd values 
 are shown in figure 7 . scatter plot of standard 
 rmsd ( x - axis ) vs symmetry - corrected 
 rmsd ( y - axis ) . each red point represents the top - scored 
 pose from systems in the sb2012 test set . the dashed lines plotted 
 by x = 2  and y = 2  
 delineate the four quadrants . the solid line indicates the diagonal . 
 the bottom insert is a magnification of quadrants i and ii . examples of 20 molecules from a total of 52 
 rescues among the top - scored 
 poses . the pdb code of the system is shown , followed by the standard 
 rmsd , and in parentheses the symmetry - corrected rmsd . each panel of figure 7 contains the 
 crystallographic 
 reference pose ( in red ) overlaid with the top - scoring docking outcome 
 ( in green ) . some of the poses , especially in the top three rows of 
 figure 7 , are easily identifiable as highly 
 symmetric molecules . the symmetry - correction algorithm has in some 
 cases improved the reported rmsd by as much as 10  , a very significant 
 change . also in figure 7 , there are examples 
 of molecules ( e.g. , 1sg0 and 1sri ) 
 that were not inverted over a major axis of symmetry , but instead 
 the rmsd was improved above the success threshold when accounting 
 for a ring flip ( as described previously ) . importantly , a breakdown 
 by protein families from the test set reveals that although the overall 
 increase in success rate is 5.0% , more dramatic success rate improvements 
 ranging up to 16.7% for specific families are observed ( table 2 ) . particularly encouraging is the 13.3% success 
 rate improvement across the large hiv protease test set ( n = 60 ) containing highly symmetric cyclic - urea based ligands , which 
 provides a challenging test case for docking . of the nine protein 
 families in table 2 for which symmetry - correction 
 appears to offer no benefit ( 0% increase in success ) , five families 
 already yield a standard rmsd success rate of greater than 50% , including 
 the highly accurate results for hiv reverse transcriptase at 95.2% 
 and neuraminidase at 93.0% . in these latter two cases , symmetry - correction the remaining four 
 families in table 2 with standard rmsd success 
 rates of less than 50% are for systems for which docking has been 
 particularly challenging ( see figure 7 in mukherjee 
 et al . ) and symmetry correction here 
 provides no improvement . as described in computational details , we have adapted the hungarian 
 algorithm to enable pruning of molecules during de novo growth to 
 avoid combinatorial explosion . the ultimate goal is a final ensemble 
 of molecules with high diversity and constructed within a reasonable 
 amount of time , that can then be used to guide organic synthesis or 
 to identify related purchasable compounds for experimental testing . 
 in terms of evaluation and validation of our de novo routines , we 
 hypothesize that hungarian - based pruning will increase efficiency 
 in experiments designed to ( 1 ) reconstruct known parent molecules 
 from restricted fragment libraries in their respective binding sites 
 ( i.e. , in less cpu time and with fewer total molecules sampled ) and 
 ( 2 ) yield increased diversity in final molecule ensembles constructed 
 from generic fragment libraries . as shown in figure 8 , the first set of experiments designed to promote 
 ligand growth was evaluated on the basis of cpu - time , number of molecules 
 processed , and tanimoto coefficient computed with respect to each 
 of 103 parent compounds . using the new hungarian algorithm - based pruning 
 heuristic , a sharp decrease in cpu time is observed ( figure 8a ) from approximately 5.4  2.6 h as the 
 number of unmatched atoms ( # -unmatched ) threshold increases from 0 
 to 10 . similarly , the total number of molecules processed shows a 
 marked decrease from approximately 15k  5k molecules over 
 the same threshold ( figure 8b ) . both of these 
 indicators demonstrate that the de novo experiments are finishing 
 more quickly and more efficiently . importantly , the average best tanimoto 
 coefficient observed during the experiment decreases only gradually 
 from 0.97  0.89 ( figure 8c ) . thus , 
 the parent molecules or analogs that are very similar to the parent 
 molecules are being rebuilt with high frequency in approximately half 
 the time , as was the goal of the experiment . ( a ) cpu time required 
 to complete the de novo experiment , ( b ) total 
 number of molecules considered during the de novo experiment , and 
 ( c ) best tanimoto coefficient to the parent compound plotted with 
 respect to the user - defined # -unmatched threshold . results are averaged 
 over all systems with 7 rotatable bonds from the sb2012 test set ( n = 103 ) . in the second set of 
 experiments , we performed de novo growth using 
 the larger seeded generic libraries ( see computational 
 details ) with a variable threshold of 05 for # -unmatched 
 ( figure 9 ) . the focus here was to test our 
 hypothesis that increasing the # -unmatched threshold in the heuristic 
 will increase diversity among the final molecule ensemble . following 
 de novo growth of new compounds , the diversity among members in each 
 ensemble was determined by computing all the possible pairwise tanimoto 
 coefficients between all molecules constructed . for example , an ensemble 
 of 50 molecules would have 50  50 pairs , minus 50 self  self - comparisons , 
 for a total of 2450 tanimoto values . for a given # -unmatched threshold , 
 tanimoto coefficients from different systems were compiled , then histogrammed 
 into discrete bins as shown in in figure 9 . histograms 
 of tanimoto coefficients for all pairwise comparisons 
 within de novo molecule ensembles . results are binned in 0.05 coefficient 
 increments over all systems with 7 rotatable bonds from the sb2012 
 test set ( n = 103 ) . figure 9 shows six final populations 
 of 
 molecules that each characterize different distributions of tanimoto 
 coefficients . as the distribution shifts to the right , molecules within 
 the ensembles are more similar to one another ( closer to 1.0 ) and 
 thus have less diversity . as the distribution shifts to the left , 
 molecules within the ensembles are less similar to one another ( closer 
 to 0.0 ) and thus have greater diversity . consistent with our expectations , 
 the peak shifted furthest to the right ( corresponding to the least 
 diversity ) occurs when the # -unmatched term is set to 0 ( figure 9 , black line ) . and , the de novo ensemble shifted 
 furthest to the left ( corresponding to the most diversity ) occurs 
 when the # -unmatched term is set to 5 ( figure 9 , purple line ) . these data indicate that , during de novo growth , 
 the hungarian algorithm is successfully identifying chemically similar 
 molecules that overlap sufficiently in cartesian space , then pruning 
 the population to make room for more diversity . in addition to the applications 
 described above , two extended applications of the hungarian algorithm 
 are envisioned : ( 1 ) as a scoring function to rank - order docked ligands 
 from virtual screening based on similarity to a known reference and 
 ( 2 ) as a tool to analyze ligand behavior during a molecular dynamics 
 ( md ) simulation . to illustrate the 
 first extended application , 
 we examined docked results from an enrichment study previously performed 
 in our laboratory , in which balius et al . docked 475 known active and 15  990 decoy compounds to the 
 target epidermal growth factor receptor ( egfr , pdb 1m17 ) . the adapted version of the hungarian algorithm for comparing 
 dissimilar molecules was used to rerank that docked ensemble by comparing 
 each docked pose to the cognate inhibitor ( i.e. , reference ) 4-anilinoquinazoline 
 ( erlotinib ) in its crystallographic binding pose with egfr . the total 
 number of unmatched atoms ( # -unmatched ) and the rmsd between the matched 
 atoms ( rmsd - matched ) were returned for each ligand . we hypothesized 
 that many of the known actives would have significant spatial and 
 chemical overlap with the erlotinib reference , and thus would be ranked 
 early in the list using the hungarian approach , thereby providing 
 good enrichment . enrichment was evaluated using four unique ranking 
 methods : ( 1 ) first by # -unmatched , and in case of a tie , then by rmsd - matched , 
 ( 2 ) first by rmsd - matched , and in case of a tie , then by # -unmatched , ( 3 ) by a combination of the # -unmatched and rmsd - matched as shown in eq 5 , and ( 4 ) the standard 
 dock grid score.5 for the ranking 
 method 
 described in eq 5 , # -ref atoms is the number 
 of non - hydrogen atoms in the reference ligand , and c1 and c2 are constants of 
 5 and 1 , respectively , chosen such that the two terms are 
 weighted approximately equally , and that greater negative values represent 
 more favorable scores . figure 10 presents the 
 receiver operating characteristic ( roc ) curves derived from each of 
 the four ranked lists along with structural overlays for top - scoring 
 compounds from two of the four methods . while visualization of the 
 roc curves can provide a qualitative assessment of enrichment , results 
 are also quantified below using area under the curve ( auc ) computed 
 across the entire database ( overall enrichment ) , and the number of 
 actives recovered after screening 1% of the database ( early enrichment ) . ( a ) 
 roc curves from docking 475 actives and 15  990 decoys 
 to egfr , then ranking the results by one of four metrics . ( b , c ) crystallographic 
 pose of reference ligand erlotinib shown as gray sticks and transparent 
 surface , identical in both panels , relative to the top four ligands 
 from the ( b ) combined value , shown as orange sticks , or the top four 
 ligands from the ( c ) dock grid score , shown as orange sticks . quantitatively , the auc data from 
 figure 10a reveal the following trend in terms 
 of overall enrichment : theoretical 
 maximum ( 1.00 ) > combined value ( 0.92 ) > # -unmatched ( 0.89 ) 
 > rmsd - matched 
 ( 0.83 ) > dock grid score ( 0.58 ) > theoretical random ( 0.50 ) . all four 
 methods provide substantially better enrichment relative to random , 
 and the three hungarian - based approaches provide enhanced enrichment 
 ( figure 10 black , red , blue lines ) relative 
 to the standard dock grid score ( green line ) , in accordance with our 
 hypothesis . in terms of early enrichment , often considered to be a 
 more useful metric for real - world applications of virtual screening , 
 the number of actives recovered at 1% of the database screened are 
 theoretical maximum ( 165 ) > combined value ( 158 ) > rmsd - matched 
 ( 131 ) 
 > dock grid score ( 113 ) > # -unmatched ( 85 ) > theoretical 
 random ( 5 ) . 
 again , all four methods provide significantly enhanced enrichment 
 relative to random , and two of the three hungarian methods provide 
 enhanced early enrichment relative to the standard dock grid score . 
 upon visual examination , the four top - ranked compounds from the combined 
 hungarian method show much better cartesian- and chemical - space overlap 
 with the erlotinib reference ( figure 10b ) compared 
 to compounds obtained using the dock grid score , which are larger 
 and have less well - overlaid groups ( figure 10c ) . overall , while the effectiveness of the scoring function should 
 be tested further on additional systems , these preliminary data indicate 
 the method shows promise in enriching for compounds that are similar 
 to a known reference molecule . as an illustration of the second extended 
 application , we explored the use of symmetry - corrected rmsd to characterize 
 ligand dynamics as a function of time . we expected that typical sampling 
 of a ligand during an md simulation of a protein  ligand complex 
 would lead to the occurrence of ring and other symmetric flips which 
 would be discernible through comparisons of standard vs symmetry - corrected 
 rmsd . figure 11 plots results derived from 
 two protein  ligand trajectories from md simulations performed 
 previously in our laboratory : ( 1 ) t4 
 lysozyme in complex with benzene ( pdb 181l ) , and ( 2 ) 
 -trypsin in complex with 1-(4-amidinophenyl)-3-(4-chlorophenyl)urea 
 ( pdb 1bju ) . here , md coordinate files saved every 1 ps were 
 least - squares fit to the starting frame ( using protein -carbons 
 as the reference ) , followed by the calculation of standard ( figure 11 , red ) and symmetry - corrected ( figure 11 , black ) rmsds for the ligand . as a practical note , 
 in the present example the ligand coordinates were extracted and then 
 postprocessed using the hungarian algorithm code in dock . in theory , 
 however , there is no barrier to incorporating the code directly into 
 an md simulation or analysis package . symmetry - corrected ( black 
 lines ) vs standard ( red lines ) heavy - atom 
 rmsds for two ligands during md simulation . in the first example ( figure 11 , left ) , 
 the benzene ligand was observed to rotate both clockwise and counterclockwise 
 about an axis normal to the plane of the ring ( as indicated by blue 
 arrows ) . although no ring  flipping  was observed in 
 this simulation , likely because of steric constraints imposed by the 
 binding site , large variations in standard rmsd were measured as a 
 direct result of the rotational motion ( figure 11 , red ) . importantly , the symmetry - corrected rmsd ( figure 11 , black ) effectively corrects for symmetry in the 
 benzene , thus highlighting variations due to either translation or 
 out - of - plane rotations , which in both instances are minimal . in the 
 second example ( figure 11 , right ) , the ligand 
 remained fairly close to its original binding configuration during 
 the md simulation . however , at approximately 550 ps , the 4-chlorophenyl 
 group flipped 180 about a rotatable bond ( as indicated by blue 
 arrows ) . corresponding to this flip , a spike in the standard rmsd 
 was observed , which was not reflected in the symmetry - corrected rmsd 
 ( figure 11 , red vs black ) . and , at approximately 
 1250 ps , the same ring flips 180 back to the original orientation , 
 and the two rmsd measures reconverge . it should be noted that although 
 this event was not easily detected visually , it was readily apparent 
 in the plotted graphs . overall , although more exhaustive tests should 
 be done , the two extended applications presented in this manuscript 
 suggest the hungarian algorithm will have use in both virtual screening 
 and ligand sampling analysis . in 
 this work , we describe implementation , validation , and application 
 of the hungarian algorithm within the dock6 program to compute symmetry - corrected 
 rmsd . the method is fast , on the order of less than a millisecond 
 per ligand ( figure 4 ) , and therefore adds only 
 a trivial amount of time to on - the - fly flexible ligand docking in 
 dock6 . more importantly , the negligible amount of time required of 
 the algorithm allows its use in de novo design . in all , 88  853 
 docked poses from the sb2012 test set of 1043 receptor  ligand 
 systems were evaluated using the new method . while the majority of 
 poses showed a < 1.0  improvement in rmsd , 
 indicative of local corrections within symmetric functional groups , 
 a smaller but significant set yielded > 3.0  improvement , 
 indicative 
 of higher symmetry corrections ( figure 5 ) . 
 a closer examination of just the top - scoring pose ( table 1 ) from each system revealed 67.7% were true positives , 
 or considered as successfully docked by either metric , and 5.0% were 
 false negatives , or considered to be rescued by the symmetry - correction 
 algorithm , for a total success rate of 72.7% . visualization 
 of specific examples of rescued poses ( figures 6 and 7 ) reveals that the 
 method is behaving as expected . highly symmetric molecules with nearly 
 perfect overlap to crystallographic references but inverted over the 
 axis of symmetry are identified and rescued , and some molecules with 
 local pockets of symmetry are also rescued . analysis of results grouped 
 by specific protein families reveals more dramatic increases in success 
 rate up to 16.7% , including the challenging hiv protease system , for 
 which an increase in success rate of 13.3% over 60 systems was observed 
 ( table 2 ) . the hungarian method was also 
 adapted to be used as a heuristic 
 for comparing spatial and chemical overlap between two dissimilar 
 molecules to aid in pruning in a de novo design version of dock undergoing 
 development . preliminary tests demonstrate that by using the heuristic , 
 sampling remains relatively unaffected , while time and efficiency 
 of computation are significantly improved ( figure 8) . in addition , diversity among molecules in the final growth 
 ensembles is higher , which for these purposes was desirable ( figure 9 ) . if preferred , it would be straightforward to 
 modify the function in order to decrease chemical 
 diversity in growing de novo populations . first , the adapted version of the algorithm was 
 used as a scoring function to rank - order compounds from a large database 
 and thus identify those related to a known reference molecule ( figure 10 ) . compounds with high overlap using the hungarian 
 metric are both spatially and chemically similar to the reference 
 compound . second , the algorithm was used to compute symmetry - corrected 
 rmsds in the course of an md simulation ( figure 11 ) . while additional testing is desirable to more fully establish 
 the utility of the algorithm for these extended applications , the accurate determination of success and failure in docking calculations 
 is essential for evaluation of new sampling routines , scoring functions , 
 protocols , and other code developments . prior to the current release , 
 the program dock contained no method to efficiently correct for symmetry 
 when computing rmsd for small molecule ligands . the source code 
 for this hungarian algorithm implementation is available to registered 
 users of dock as part of the latest release at http://dock.compbio.ucsf.edu .
false 
 negative docking outcomes for highly symmetric molecules 
 are a barrier to the accurate evaluation of docking programs , scoring 
 functions , and protocols . 
 this work describes an implementation of 
 a symmetry - corrected root - mean - square deviation ( rmsd ) method into 
 the program dock based on the hungarian algorithm for solving the 
 minimum assignment problem , which dynamically assigns atom correspondence 
 in molecules with symmetry . 
 the algorithm adds only a trivial amount 
 of computation time to the rmsd calculations and is shown to increase 
 the reported overall docking success rate by approximately 5% when 
 tested over 1043 receptor  ligand systems . for some families 
 of protein systems the results are even more dramatic , with success 
 rate increases up to 16.7% . 
 several additional applications of the 
 method are also presented including as a pairwise similarity metric 
 to compare molecules during de novo design , as a scoring function 
 to rank - order virtual screening results , and for the analysis of trajectories 
 from molecular dynamics simulation . 
 the new method , including source 
 code , is available to registered users of dock6 ( http://dock.compbio.ucsf.edu ) .