in drug discovery projects , the three - dimensional ( 3d ) structure of the receptor is not always available . in such cases , pharmacophore model building , 3d quantitative structureactivity relationship ( qsar ) , comparative molecular field analysis ( comfa),(1 ) and ligand - based virtual screening all depend on a good algorithm to flexibly align small molecules . in light of the importance of the subject on top of the methods mentioned in the comprehensive review of lemmen and lengauer,(2 ) new efforts have continued to appear . similar to docking , the small molecule alignment problem can be divided into two parts . second , a scoring function needs to be able to distinguish good alignment poses from other possibilities . the current work is about obtaining a good scoring function for this purpose . in principle , a scoring function that can distinguish the correct alignment pose from incorrect ones is not necessarily the same as a scoring function that upon minimization gives poses that resemble the correct one as much as possible . in practice , a good scoring function probably performs well for both purposes . in this work , we are primarily interested in a scoring function for distinguishing good poses . in any case , we believe that the exact answer to the molecular alignment problem may not be simple . for example , using only  carbons in the binding pocket for superposition may give a slightly different answer from that using all the pocket atoms or using the pharmacophore elements involving the ligand and the pocket . moreover , the ligand may have some freedom of movement inside a pocket . hence trying to nudge down the geometric difference between a proposed alignment hypothesis and the  correct  alignment all the way to zero may not be necessary , or even meaningful , as long as the correct alignment mode is obtained . for scoring functions used in docking , one can use molecular mechanics force fields derived from first principles . such functions need to be derived from some statistical analysis of a training set of known molecular alignments . the protein data bank(17 ) ( pdb ) is an obvious source to obtain such a training set . just in the five years since 2003 , the number of entries has more than doubled ( see http://www.pdb.org ) . moreover , many old entries have been revisited and cleaned up , and many low - resolution structures have been superseded by better ones . by now , there are hundreds of systems of entries with identical proteins but different ligands . we believe that now is a good time to make use of the growth in the pdb to revisit the derivation of a scoring function for molecular alignment . for the problem of small molecule alignment , there are two main types of scoring functions . the first type is atom based : when two molecules are being aligned , the score consists of a sum of terms that are based on intermolecular atom pairs ( i.e. , each pair has one atom coming from each molecule ) . the second type of score is field based : the electrostatic or steric fields of the molecules or their surfaces are compared to arrive at a score . this would be somewhat slower to compute since the fields or the molecular surfaces need to be calculated from the atomic coordinates . hence it would be difficult to flexibly refine an ensemble of aligned molecules or to deal with a huge number of conformers . indeed several of the field - based alignment methods rely on an independent conformational search engine to generate good conformers ( e.g. , shapelets,(6 ) brutus,(13 ) the  molecular field extrema  method,(11 ) and mimic(18 ) ) . since our scoring function will be used to distinguish good binding poses among a huge number of possibilities , it needs to be calculated quickly . an atom - based function would hence be more appealing than a molecular field - based function . when creating a hypothesis for an alignment of molecules , it is desirable for atoms ( from different molecules that are being aligned ) with similar properties to be near to each other . intermolecular atom pairs that are far away should have little influence on the goodness of an alignment . this can be achieved by associating each intermolecular atom pair with similar properties with a gaussian scoring term : w exp(r ) , where r is the distance between the two atoms ,  is a parameter controlling the range of the interaction , and w is the weight of this term . it is also desirable for the alignment hypothesis for atoms with incompatible properties ( e.g. , hydrophobic vs hydrophilic ) to not superpose . this can be achieved using the same functional form but with a negative value for w. besides atoms , hydrogen - bond projected features can also be considered . these features represent the expected positions of the hydrogen - bonding partner on the binding pocket based on the ligand geometry . the overlap of hydrogen - bond projected features of the same type should be favorable for an alignment hypothesis . again , the same gaussian score can work , with r being the distance between the intermolecular features . while the gaussian function is widely used in molecular alignment methods ( e.g. , fap,(4 ) pharao,(5 ) moe,(16 ) sq,(19 ) and seal(20 ) ) , there are other functional forms that can provide a similar effect . these include piecewise linear functions ( e.g. , flame),(12 ) and piecewise linear and quadratic functions ( e.g. , tfit).(21 ) it can be construed that indeed any bell - shaped function can work in a similar fashion . in their pioneering work 20 years ago , kearsley and smith(20 ) noticed that the lorentzian function ( 1/(1 + r ) ) yielded similar results to the gaussian but was significantly faster to compute . however , the advance in computational technology to calculate exponents has eliminated this advantage of the lorentzian function . for this work we chose the gaussian because it involves a minimal number of tunable parameters and has a simple , continuous shape . another consideration is that function optimization procedures often require the first and second derivatives of the function . in the case of the gaussian , once the function has been calculated , its first and second derivatives can be obtained without much extra work . a scoring function s could thus be constructed based on a sum of gaussian terms : 
  in the formula , tk is a term encouraging / discouraging the overlap of a particular pair of types ( ak , bk ) of atoms / features ; wk is the weight of the term ; i and j are atoms / features from the two molecules ; rij is the distance between the intermolecular atom / feature pair ; and  is a parameter controlling the range of this type of interaction . ligands binding in the same mode are expected to have substantial volume overlap . the first term in the scoring function , t0 the terms t1 and t2 encourage the overlap of pairs of intermolecular donors and acceptors . after these  attractive  terms come the  repulsive  terms that discourage the overlap of intermolecular atom pairs of incompatible types . first , hydrophobic and hydrophilic atoms ( atoms capable of hydrogen bonding ) are not expected to occupy the same part of the binding pocket , as reflected by t4 and t5 . second , hydrogen - bond acceptors that are not also donors are expected to avoid positions of donors that are not also acceptors . the last two terms , t8 and t9 , are based on hydrogen - bond projected features . the scoring function consists of a summation of terms , tk s , each of which has the form : 

 each row of the table corresponds to one term , tk , and i and j denote atoms ( for t0t7 ) or features ( for t8 and t9 ) coming from the first and second molecule , respectively . these atoms / features have to be of certain types , ak and bk , as given in the second and third columns of the table , and rij is the distance between the atoms / features . the names of the various constants for the different terms , wk and  , are given in the fourth and fifth columns of the table . all the atom - based terms have the same value for  , namely a , while all hydrogen - bond projected feature - based terms have another value for  , namely p . to keep the number of parameters low , the same value of  , a , the projected feature terms used a different value of  , p , allowing these terms to be more diffuse than the atom - based terms if necessary . atom - based terms were followed by projected feature terms . for the atom - based terms , based on this order , the scoring function was built up term by term . at each stage , the one or two parameter(s ) of the newly added term was / were optimized . a series of scoring functions were created using different values for the newly added parameter(s ) . the parameter value(s ) corresponding to the best scoring function in the series would be picked . after the scoring function was optimized , it was incorporated into the flexible alignment functionality in molecular operating environment ( moe).(22 ) molecular alignment results were compared with two recent publications . proschak et al.(6 ) performed rigid - body alignments for eight thermolysin ligands using their shapelets method and compared their result with two previous alignment methods . we performed validation runs on these thermolysin ligands to see how our scoring function fared . chen et al.(23 ) studied the geometric accuracy of aligning small molecules using rocs(26 ) and flexs.(27 ) they used eight protein targets as test sets : hiv protease , cyclin - dependent kinase 2 ( cdk2 ) , estrogen receptor 1 ( esr1 ) , mitogen - activated protein kinase 14 ( p38 ) , thermolysin , human rhinovirus capsid , elastase , and trypsin . six of these ( all but thermolysin and elastase ) consist exclusively of structures from the publicly accessible pdb.(17 ) we used these six systems for validation runs . both rigid - body alignments and flexible alignments were compared . for training a scoring function for the validation runs , all systems that contain any pdb entries involved in the test sets were naturally excluded . ligands from identical or almost identical pockets ( see below for details ) from the pdb were used for training the scoring function . only pdb entries with a resolution of better than 2.5  were used . ligands were defined as freestanding molecules with no more than 300 heavy atoms and 10 residues . in order to exclude waters and very small molecules , each ligand must have at least one pair of heavy atoms separated by three bonds . ligands without carbon atoms or with exotic elements ( elements other than c , h , n , o , s , p , f , cl , br , and i ) were discarded . atp analogs were excluded because their large population in the pdb would likely skew the data . ligands without any rings were also excluded so as to remove the lipids without excluding any drug - like compounds . we first obtained a list of protein domains from structural classification of proteins ( scop).(28 ) for each domain , all pdb entries containing that domain were superposed using the protein backbone of that domain . overlapping ligands were clustered into groups so that each group would correspond to one pocket . for each pocket , since we started from scop domains , it was possible that pockets within a group could involve different residues or have different geometries . for our purpose of scoring function training , it may not be very meaningful to study the alignment of ligands in different pockets . a table was then constructed with each row corresponding to one protein and each column corresponding to one  carbon cluster . a cell of the table has the residue type as the entry if the protein contributes to that cluster . the pocket similarity between a pair of proteins was defined to be the fraction of table columns that are identical ( i.e. , same residue type , or empty cell , for all cells in the column ) . after the pocket similarities between all protein pairs were determined , the proteins were clustered , using complete linkage and pocket similarity as the metric , with a cutoff of 0.8 . in other words , any two pockets in the same group had at least 80% of the table columns identical . each original pocket group might thus be divided into multiple groups . for each resulting group , a fine superposition on the 3d structures was performed again using  carbons common to all pockets , since the definition of the pocket might have changed due to the change in the content of the ligand set . this whole procedure was repeated until there was no new division of the pocket groups . upon inspection , some systems were found to be uninteresting in the sense that the aligned ligands only displayed one type of structure at most regions of the pocket . studying this type of system would not give us more information than studying self - alignments . . therefore the training ligands were restricted to those having no more than 12 rotatable bonds . finally , to ensure diversity , a set of representative ligands was chosen for each system such that no pair of representatives had a tanimoto similarity of over 0.8 on their maccs fingerprints . only 100 systems remained after these final filters , and 85 systems had no more than 20 representative ligands . this is our comprehensive aligned small molecules data set . for training a scoring function for the validation runs , 11 of the 100 systems in the comprehensive data set were excluded because they contained entries that would be used in the runs . this was more than the number of test systems ( 7 ) because some test systems ( e.g. , the 57 ligands in the cdk2 system ) had variations in the pocket residue composition or geometry and were hence subdivided into more than one pocket group when our comprehensive data set was constructed . to keep things simple , the training of the scoring function was based on trying to align a flexible source ligand onto a rigid target ligand . dummy poses were generated by making small adaptations to the moe docking engine,(22 ) as described below . n aligned ligands , all n  n pairs of ligands were considered with the following requirement . the fraction of heavy atoms in the target that lie within 2  of any heavy atom of the source must exceed 0.5 , and the fraction of heavy atoms in the source that lie within 2  of any heavy atom of the target must exceed 0.75 . this excluded ligands with different binding modes and ensured that no source ligand was substantially larger than the target . if this resulted in more than 5000 conformers , a random subset of the systematically generated conformers would be used . only conformers with severe atomic overlaps this is because we envisioned that when the scoring function is used to distinguish poses with the correct binding mode , any poses that are fished out can be refined to improve the pose geometry and any atomic overlap within the molecule will be relieved . for each sourcetarget pair , each time , a random set of three heavy atoms on the target ligand was chosen , together with a random conformer of the source ligand . triplets of heavy atoms on the source ligand , which formed a triangle similar to that formed by the three chosen target atoms , were used to align the source ligand onto the target ligand . poses whose root - mean - square deviation ( rmsd ) for the heavy atoms was under 1.5  were considered duplicated , and only one copy was kept . for each scoring function , the scores for all dummy poses were calculated . for each set of alignments involving a particular pair of source and target ligands , we counted the number of  good  poses , n25 , among the top 25 scoring poses . as for whether a pose is  good  or not , we used fuzzy logic . poses with an rmsd ( heavy atoms only ) from the correct pose of under 2.0  were considered good . poses with an rmsd between 2.0 and 3.0  were considered partial good poses , with the  partiality  decreasing linearly from 1 to 0 . for each system ( i.e. , same pocket ) , we averaged the number of good poses for all sourcetarget pairs in the system . the receiver operating characteristic ( roc ) curve is often used to measure the success of scoring functions for virtual screening.(29 ) in short , the roc curve plots the fraction of true positives recovered ( out of all active compounds ) against the fraction of false positives obtained ( out of all inactive compounds ) , as the cutoff value varies . it is generally considered that the bigger the area under the curve , the better is the scoring function . we felt that , for our purpose , using the area under the roc curve was not as suitable as using the number of good poses amongst the top 25  , since we generally can only afford to consider a limited number of top scoring poses . hence all but the very beginning part of the roc curve would be of significance to us . we think that human inspection can typically handle 10 or so poses . on the other hand , if a thorough refinement is desired as part of an automatic process , then it would be reasonable to intensively process up to about 100 poses per ligand . hence here we used the  top 25 poses  . as a precaution , we repeated all the calculations and validation runs using the  top 100 poses  instead of the  top 25  . it turned out that only one out of the six optimized parameters for the scoring function would come out different , and the quality of the of validation runs remained very similar . the details are given in section a of the supporting information . as mentioned previously , the scoring function was built up term by term . at each stage , only one ( or two ) parameter was introduced . each table corresponds to one stage of optimization where one ( or two ) parameter was optimized . the resulting numbers of good poses among the top 25 , n25 , corresponding to the various values of the parameter being optimized are given in the tables . the optimal value of the parameter and its corresponding value of n25 are highlighted in bold . this table gives the average number of good poses , n25 , within the top 25 scoring poses for various values of a . we started by including only the volume overlap term ( t0 in table 1 ) in the scoring function . in other words , to a first approximation , the bigger the volume overlap between two molecules , the better was the score . in order to make the scoring function quick to calculate , the only parameter to be optimized was . obviously a smaller value for  would result in a more diffuse gaussian function . consider the simple case to align one atom ( the source ) onto a diatomic molecule ( the target ) . this is similar to trying to align methane , with one heavy atom , onto ethane , with two heavy atoms . suppose the internuclear distance of the diatomic molecule is 1.5  ( for comparison , a cc bond in a benzene ring is about 1.40  , and a cc bond in an ethane is about 1.54  ) . the shapes of the scoring function for three different values of  are plotted in figure 2 . for large values of  , the source atom will see a score with two distinct peaks located very close to the two target atoms . it will have a strong tendency to align with either one of the target atoms . for example , if  is 2.0  ( the blue curve of figure 2 ) , then the distance between the peaks will be 1.46 . the value of the score at the trough midway between the peaks will be 64% of the height of the peaks . when  decreases to 1.2  ( the black curve of figure 2 ) , the peaks will move toward each other , away from the atomic positions . the distance between the peaks will drop to 1.18  , and the trough to peak height ratio of the score will increase to 94% . in other words , the score has already very much smoothed out the atomic positions . the source atom can settle for any position between the two atoms with little change on the score . how the scoring function is controlled by the parameter  , for the case of a target molecule consisting of two atoms , located at positions 0.75 and 0.75 . for a source atom to be superimposed onto the target molecule , the scoring function is exp(r12 ) + exp(r22 ) , where r1 and r2 are the distances from the source atom to the two target atoms . the scoring functions corresponding to  values of 2.0 , 1.5 , and 1.2 are plotted in blue , red , and black , respectively . given a scoring scheme and an initial aligned pose between two rigid ligands , one can refine the alignment by optimizing the score . when we did this with many starting poses , we found that we ended up with fewer resulting aligned poses when we used smaller values of  ( unpublished results ) . this is in line with our expectation that smaller values of  smooth out the score landscape . according to previous works of gaussian function - based alignment scores , values of  of between 0.1 and 0.5 , we tried out a total of 7 values for a , from 0.06 to 0.7  . the corresponding numbers of good poses among the top 25 are tabulated in table 2a . it can be seen that for the values of a tested , a value of 0.5  was the best , with an average of 4.83 good poses among the top 25 . after establishing a value for a , the relative weights of the atom - based terms were optimized one by one . we first looked at the weight wda for the hydrogen - bonding atoms ( donor and acceptor ) similarity terms ( t1 and t2 in table 1 ) . for this we tried out 7 values for wda , from 0.5 to 32 , in logarithmic intervals . the corresponding numbers of good poses among the top 25 are tabulated in table 2b . a value of 4 turned out to be the most favorable , with an average of 6.10 good poses among the top 25 . the new terms are t1 and t2 , which correspond to the attraction between hydrogen - bond donor / acceptor atoms . this table gives the average number of good poses , n25 , within the top 25 scoring poses for various values of wda . after this , we turned to the similarity term between hydrophobic atoms ( t3 in table 1 ) . besides trying to establish the optimal weight whh for this term , we also tried out two schemes of defining hydrophobic atoms . one used the conventional hydrophobic atom assignment in moe.(22 ) the other scheme was based on how far topologically an atom was from any hydrogen - bonding atom . atoms that were two or more bonds away from any hydrogen - bond acceptors or donors were considered hydrophobic . tables 2c and 2d give the numbers of good poses among the top 25 for various values of whh using the two schemes . it can be seen that the two schemes gave very similar results , with the optimal weighting factor whh to be around 1.0 for both schemes . the scoring function at this stage is ( t0 + t1 + t2 + t3 ) . this table gives the average number of good poses , n25 , within the top 25 scoring poses for various values of whh . the scoring function at this stage is ( t0 + t1 + t2 + t3 ) . this table gives the average number of good poses , n25 , within the top 25 scoring poses for various values of whh . next we considered the repulsion term between hydrogen - bonding atoms and hydrophobic atoms ( t4 and t5 in table 1 ) . here we see that the first scheme of defining hydrophobicity was slightly better and that the optimal weight for the repulsion term , wr1 , was around 1.0 . the scoring function at this stage is ( t0 + t1 + t2 + t3 + t4 + t5 ) . the new terms t4 and t5 correspond to the repulsion between hydrophobic and hydrogen - bonding atoms . this table gives the average number of good poses , n25 , within the top 25 scoring poses for various values of wr1 . the scoring function at this stage is ( t0 + t1 + t2 + t3 + t4 + t5 ) . the new terms t4 and t5 correspond to the repulsion between hydrophobic and hydrogen - bonding atoms . this table gives the average number of good poses , n25 , within the top 25 scoring poses for various values of wr1 . after this we investigated a repulsion term between hydrogen - bond donor and acceptor atoms ( t6 and t7 in table 1 ) . perhaps surprisingly , it was found that this term did not seem to improve the scoring function ( see table 2 g ) . the scoring function at this stage is ( t0 + t1 + t2 + t3 + t4 + t5 + t6 + t7 ) . the new terms t6 and t7 correspond to the repulsion between hydrogen - bond donor atoms that are not acceptors and hydrogen - bond acceptor atoms that are not donors . refer to table 1 for the details of these terms . this table gives the average number of good poses , n25 , within the top 25 scoring poses for various values of wr2 . finally the projected acceptor and donor feature points were considered ( t8 and t9 in table 1 ) . these points correspond to the expected positions of donors and acceptors on the binding pocket . in this sense so we tried out different values of  as well as different weights , wpf , for this term . it can be seen that an p value of 0.125  coupled with a weight of 2.0 was a good combination . otherwise an p value of 0.0625  coupled with a weight of 1.0 was also respectable . the scoring function at this stage is ( t0 + t1 + t2 + t3 + t4 + t5 + t8 + t9 ) . the new terms t8 and t9 correspond to the attraction between hydrogen - bond donor / acceptor projected features . refer to table 1 for the details of these terms . this table gives the average number of good poses , n25 , within the top 25 scoring poses for various values of wpf and p . to summarize , it was found that a good scoring function is of the form given by the terms of table 1 , except terms t6 and t7 , with the values of the various coefficients given in table 3 . the optimal scoring function consists of the terms t0t5 , t8 , and t9 , as given in table 1 , with the various parameters as given by this table . for each of the parameters considered however , there was a considerable number of ligandpocket systems , and in each system , there was a considerable number ( at least 25 , but often much more ) of pairs of sourcetarget ligands . more importantly , the fact that there are smooth trends in the numbers in tables 2a2h , as the parameters were varied , serves to confirm that the numbers are accurate enough for our purpose of locating the optimal values for the parameters . moe s flexible alignment ( flexalign ) functionality(22 ) was used for the validation runs . for rigid - body alignments , conformers from the pdb flexalign generated random poses of the source molecule superposed onto the target molecule and optimized the alignment score s , while keeping the molecules rigid . for the flexible alignments , for each source molecule , one randomized conformation was used as input . flexalign generated a random conformer from this input , randomly superposed it onto the target molecule , and then optimized this pose under the influence of the alignment score between the source and target molecules . the target molecule remained fixed throughout the process . in order to avoid unrealistic conformations , the internal energy u of the source molecule must also be considered alongside the alignment score s. for the internal energy , we used the mmff94 force field with a distance - dependent dielectric . pose optimization was carried out by minimizing the effective energy , e , which was defined by 
  alignment solutions were also ranked according to their e values . a value of 4 kcal / mol was used for k. this value seemed to give a good balance between the internal energy u and the alignment score s. the difference in the force field energy u between the conformation of the source molecule in the top scoring ( i.e. , the one with the lowest e value ) solution and its corresponding local energy minimum was generally small . the average of this difference ( over the n  n cross - alignments for each system ) for the cdk2 , hiv , p38 , esr1 , trypsin , and rhinovirus systems were , respectively , 0.32 , 0.76 , 0.25 , 0.33 , 0.25 , 0.25 kcal / mol per heavy atom . in contrast , the average of the difference in the internal energy u between the ligand conformer in the crystal structure and its corresponding local energy minimum is over 1 kcal / mol per heavy atom for all the six systems . we also investigated whether the value of the alignment score s or the effective energy e gives an indication of the goodness of the alignment . our results were first compared with those in table 4 of proschak et al.(6 ) this involved pairwise rigid - body alignments on a system of eight thermolysin ligands . for each ligand pair , one ligand was used as the target , while the other was used as the source . an alignment was considered successful if the top scoring pose had a heavy atom rmsd of under 2  from the correct answer . our results , together with those obtained by proshak et al.,(6 ) are tabulated in table 4 . according to proshak when surfcomp(24 ) was used with the electrostatic potential , the success rate was 75% . in contrast , our scoring function yielded a success rate of 93% ; there were only 2 failed cases out of 28 alignments . our scoring function returned the correct answer as the second best scoring solution ( figure 3 , bottom ) . in any case , our best scoring solution looks highly plausible with heavy overall volumetric overlap as well as the overlap of the aromatic ring and two carbonyl groups ( figure 3 , top ) . carbons of our solutions ( top and bottom are best and second best scoring , respectively ) are in cyan . the eight thermolysin ligands have pdb i d s 1thl , 1tlp , 1tmn , 3tmn , 4tmn , 5tln , 5tmn , and 6tmn . for each pair of nonidentical ligands , one ligand was chosen as the template , and the other ligand was aligned onto the template . if the top scoring solution yielded an rmsd of under 2  with the crystallographic pose for the heavy atoms , then it was considered successful . this table gives the success rates for the 28 cases , as reported on table 4 of proschak et al.(6 ) tests were also run on six systems mentioned in chen et al.,(23 ) where rocs(26 ) and flexs(27 ) were used . an alignment was considered correct if the top scoring solution had a heavy atom rmsd of under 2  with the crystallographic overlay . except for the rhinovirus system , which will be examined in more details below , and for the flexible alignment of the cdk2 system , our results are significantly better than those obtained by chen et al . using flexs or rocs . for rigid - body alignments , our fraction of correct results was at least 40% and reached 80% or more in two systems . naturally , flexible alignment is more difficult than rigid - body alignment , and so the percentages of correct results are lower for all systems except rhinovirus ( see analysis below ) . for each of the six systems , pairwise cross - alignments were carried out for all n  n ligand pairs , where n is the number of ligands in the system . if the top scoring result yielded an rmsd of under 2  for the heavy atoms , then it was considered correct . this table gives the percentage of correct results for the n  n ligand pairs in the system , as reported on table 2a of chen et al.(23 ) for the cdk2 system , our results for the rigid - body alignments but not for the flexible alignments are significantly better than those of chen et al . cdk2 is a large system ( 57 ligands ) with a complex pharmacophore model.(31 ) if the six ligands 1di8 , 1p5e , 1pkd , 1pxi , 1pxj , and 1wcc were excluded , then our results would improve from 39.6% to 47.6% correct for the rigid - body alignments and from 22.4% to 26.6% correct for the flexible alignments . of the 57 ligands in the system , only 7 ligands have more than 33 or less than 17 heavy atoms , and 1p5e , 1pkd , 1pxi , 1pxj , and 1wcc are among them , with 13 , 36 , 14 , 14 , and 8 heavy atoms , respectively . since the pharmacophore model is complex,(31 ) it seems that volumetric overlap is an important factor for alignment . hence ligands that are much smaller or much larger than average would be difficult to align . moreover , an important pharmacophore element is the set of two hydrogen - bond donors and one acceptor interacting with the hinge region of the protein . these correspond to pharmacophore features d1 , d2 , and a1 mentioned in zou et al.(31 ) and are shown in figure 4 . out of the 57 ligands , only 4 do not have at least 2 out of these 3 features , and ligands 1di8 , 1p5e , and 1wcc are among them . important pharmacophore features for the cdk2 ligands include two hydrogen - bond donors ( magenta sphere ) and one acceptor ( cyan sphere ) . the fact that they are relatively comparable in size ( all but three have between 39 and 61 atoms ) , do not have a homogeneous ( spherical ) shape , and have rich feature sets might have helped the rigid - body alignment result , because there are not that many possibilities for the molecules to align well . one possibility for our superior results is that we have trained our scoring function using a significantly larger data set from the pdb than the version of flexs or rocs that were used . of all six test systems , the drop in performance in going from the rigid - body to the flexible mode is most pronounced for the hiv system . of the 28 hiv ligands , 5 have 10 rotatable bonds . 7 have between 12 and 16 rotatable bonds . for the 13 ligands of the p38 system , our scoring function produced an overall correct rate of 43% and 30% for the rigid - body and flexible alignments , respectively ( table 5 ) . it is worth noting that one of the ligands , 1wbo , is significantly smaller than the others ( figure 5 ) . it can be seen that the ligand 1di9 ( magenta in figure 6 ) binds in a mode different from the others . the four ligands 1m7q , 1ouk , 1ouy , and 1ove ( yellow in figure 6 ) also bind in a mode somewhat different from the others . for the remaining 7 ligands , our scoring function actually got 73% of the pairwise cross - alignments correct in the rigid - body mode , and 53% correct in the flexible mode . and for the 16 pairwise cross - alignments of the four ligands 1m7q , 1ouk , 1ouy , and 1ove , our scoring function got them all , 100% , correct in the rigid - body mode and 94% ( all but one of the 16 ) correct in the flexible mode . ligands within each group have comparable sizes and can be considered to bind in a similar mode , as can be seen in figure 7 . for the pairwise cross - alignments of all 13 ligands , our scoring function got 59% and 41% correct in the rigid - body and flexible modes , respectively . however , for the group of 6 smaller ligands , we got 94% of the 36 rigid - body cross - alignments correct and 72% of the flexible alignments correct . for the group of 7 larger ligands , we got 92% of the 49 rigid - body cross - alignments correct and 65% of the flexible alignments correct . the 6 smaller ligands ( 1a52 , 1gwq , 1l2i , 1x7e , 1x7r , and 3erd ) are in cyan . the 7 larger ligands ( 1r5k , 1sj0 , 1uom , 1xp1 , 1xp9 , 1xqc , and 2bj4 ) are in magenta . the crystallographic alignment of the 7 trypsin ligands is given in figure 8 . ligand 1tni ( cyan in figure 8) assumes a binding mode different from the others . for the remaining 5 ligands , our scoring function got all 25 cross - alignments correct in both the flexible and rigid - body modes . figure 9 gives our result for flexibly aligning 1tnj onto the much larger ligand 1pph . our second best scoring solution ( figure 9 , bottom ) yielded the correct alignment mode . as for our best scoring solution ( figure 9 , top ) , it involves a perfect atom to atom matching for all atoms of 1tnj . the aromatic ring is mapped to an aromatic ring of 1pph , all aliphatic carbons are mapped to aliphatic carbons on 1pph , and the nitrogen atom is mapped to a nitrogen on 1pph . our third best scoring solution ( figure 10 , bottom ) roughly yielded the correct binding mode . our best scoring solution ( figure 10 , top ) involves a matching of atoms similar to that of the best scoring solution of aligning 1tnj onto 1pph ( figure 9 , top ) . our second best scoring solution ( figure 10 , middle ) also involves a good matching of atoms , with the phenyl ring of 1tnl aligned to a six - membered ring of 1pph . the remaining ligands are 1tnh , 1tnj , 1tnk , 1tnl , 3ptb . flexibly aligning 1tnj onto 1pph ( carbons in green ) . carbons of our solutions ( top and bottom : best and second best scoring ) are in cyan . flexibly aligning 1tnl onto 1pph ( carbons in green ) . carbons of our solutions ( from top to bottom : first , second , and third best scoring ) are in cyan . for the 8 rhinovirus ligands , our scoring function got 50% of the pairwise cross - alignments correct in both the rigid - body and flexible modes . as shown in figure 11 , all eight ligands have a long shape and are almost symmetric . there is a heterocyclic ring on either end of each ligand , connected by a long , flexible linker chain . as noted in the original publication for the x - ray structures,(32 ) the ligands have two binding modes , one being the inverse of the other , as shown in figure 11 . it is easy for the alignment algorithm to superpose the two rings on either end but difficult to distinguish between the two binding modes . hence all the cross - alignments within a binding mode were obtained correctly , but the inverse binding mode was obtained for each pair across the different binding modes . the carbons are in magenta for one mode ( 2rm2 , 2rr1 , 2rs1 , and 2rs3 ) and in cyan for the other mode ( 2r04 , 2r06 , 2r07 , and 2rs5 ) . when the scoring function was trained for the validation runs , systems involved in the test sets were excluded . after the validation runs , the training process was repeated to make use of all data , including systems involved in the test sets . it turned out that using all the data did not lead to any change in the scoring function . the updated tables 2a2h are given in tables ah of section b of the supporting information . although the numbers of good poses among the top 25 ( n25 ) have changed slightly , the positions for the optimal values for all the parameters have not . based on the pdb(17 ) and scop,(28 ) a comprehensive set of aligned ligands binding in the same pocket has been compiled . we used only high - resolution structures and considered only pockets that are highly similar and have diversified , drug - like ligands . for each pair of ligands binding in the same pocket , one ligand was used as the target , and dummy aligned poses of the other ligand were generated . by studying the correlation between the score of these dummy poses and their similarities with the correct answer , a small molecule alignment scoring function was built up term by term . eventually we obtained the scoring function composed of terms given by table 1 , except terms t6 and t7 , with parameters as given by table 3 . we verified our scoring function by comparing results of rigid - body and flexible alignments performed using the flexible alignment functionality in moe(22 ) with those reported in the literature . we found that our results are superior to those of five methods reported in two recent publications . we are happy to share our comprehensive aligned small molecules data set with other scientists so that they can use it to optimize their own scoring functions .