BACKGROUND
drug design and discovery play a pivotal role in driving research in computational chemistry and biology  <cit> .
in computational drug design and discovery, it is often necessary to determine, as a first step, the binding of a ligand to a target protein.
the computational scheme for predicting ligand binding occurrence, affinity, and orientation is commonly referred to as "molecular docking", which has been a topic of intensive research for decades  <cit> .
the development of a molecular docking tool usually starts with an efficient search algorithm, which places the ligand in the active site of the target protein in numerous different positions, orientations, and, in flexible docking, conformations.
these are then evaluated by a scoring function to distinguish between good  and bad  docking solutions.
the two aspects of searching and scoring can be, and usually have been, developed and evaluated separately, although one clearly affects the other and a balance is often sought to meet specific study aims  <cit> .
scoring functions for molecular docking are traditionally either physics-based or knowledge-based  <cit>  and differ mainly in the derivation of the mathematical models used to compute the energies of molecular interactions.
physics-based scoring functions employ a set of molecular interaction terms to compute binding energies.
for example, the scoring function used in g-score  <cit>  or autodock  <dig>   <cit>  is based on the molecular mechanics force field used in tripos  <cit>  or amber  <cit> , respectively, and f-score  <cit>  and chemscore  <cit>  derive the coefficients of their energy terms using regression analysis of experimentally determined binding energies  <cit> .
in contrast, knowledge-based scoring functions, such as pmf  <cit>  and drugscore  <cit> , rely on statistical observations of preferred protein-ligand contacts, from which binding energies are calculated.
a common thread in these scoring functions is that, with few exceptions, they all operate on the assumption that the total interaction will be faithfully represented by an additive summation over a series of pairwise interactions between interaction centers, which are either real  or imaginary .
despite fundamental fallacies in this assumption, it has been adopted for many years, as otherwise the problem would be intractable or too computationally expensive to handle practically  <cit> .
in this study, we investigated a new idea for developing a scoring function for molecular docking.
the novelty of the idea lies in the use of a network approach to identify frequently occurring patterns of protein-ligand interactions not in the form of pairwise interactions, but in the form of network motifs, where "network" refers to a collection of pairwise interactions.
the motifs, which are assemblies of multiple pairwise interactions between proteins and ligands that are frequently observed in a database of known protein-ligand complex structures, thus represent energetically favorable ways of positioning a ligand molecule in the active site of a protein.
such a motif-based approach offers a new line of exploration for research on molecular docking.
as a first step in this line of research, we present the construction of a motif-based scoring function, called motifscore, and the evaluation of its ability to distinguish between good and bad docking solutions.
our results using a benchmark dataset showed that motifscore performed well against a number of scoring functions used in existing docking programs.
methods
in constructing the protein-ligand interaction network, the three-dimensional  coordinates of a set of protein-ligand complexes were each transformed into an atom-atom interaction network.
there are two types of nodes in the network, namely atoms  from the protein and the ligand.
two nodes  are connected by an edge if their interaction is deemed significant by a distance threshold determined from a statistical analysis of a dataset of protein-ligand complexes .
interactions between atoms within the protein or within the ligand are not considered.
these networks were then broken down into many interaction motifs representing simple units of the network with specific protein-ligand interaction patterns.
datasets
to establish distance thresholds for constructing the protein-ligand interaction network, we performed a large-scale statistical analysis on a diverse set of known structures of protein-ligand complexes.
the available protein-ligand complex structures as of january  <dig>   <dig> in the protein data bank   <cit>  were screened.
to simplify our task, only those determined by x-ray crystallography and with only one ligand molecule were selected.
furthermore, those containing dna/rna molecules and those in which the ligand binds to its target protein covalently were discarded.
complexes containing a heme group were also excluded because one may consider the heme group as a part of the protein, rather than a ligand.
solvent molecules and ions, such as chloride or ammonium, that are often included to facilitate crystallographic studies were not considered as ligands, whereas metal ions such as zinc or calcium ion that are coordinated by metal-binding amino acids were considered as an extension of the protein molecule, and their interactions with ligand were considered part of the protein-ligand interaction network constructed.
in all, we obtained a total of  <dig>  structures of complexes that could be used in the statistical analysis.
this set of pdb protein-ligand complex structures is not entirely non-redundant since a single protein or its mutants may be co-crystallized with many similar compounds, and, conversely, the same compound may be co-crystallized with multiple homologous proteins.
however, unless both the protein and the ligand are completely identical, part of the resulting interaction networks can still be distinct, and as defining non-redundant networks is not a straightforward task, we chose to include as many network connectivities  as possible while conducting normalizations to minimize the potential bias of using a redundant dataset.
furthermore, analysis on a reduced set  by excluding those that share >25% protein sequence identity and identical ligand showed that the resulting network motifs and their grades  had a pearson's correlation coefficient of  <dig>  with those of the set of  <dig>  structures , suggesting that including homologous proteins and identical ligand compounds would not significantly affect the outcome of our analysis.
besides the set of  <dig>  structures, we used a subset of the lpdb   <cit> , which contains decoys generated from molecular dynamics simulations  <cit> , to optimize the parameters of our scoring function, motifscore.
there are more than  <dig> complexes, each with a set of various numbers of decoys, in the lpdb, but only  <dig> of these were used in this work because we found that the ligand names, atom names, and the ordering of atoms of many of the decoys in lpdb were not the same as in the original data files in the pdb, making it difficult for us to correctly assign their designated atom types .
for testing and evaluating motifscore, the data set of wang et al.
<cit>  was used.
the wang dataset consists of a set of  <dig> protein-ligand complexes, each of which comes with  <dig> docked conformations  generated using autodock  <dig>   <cit> .
since wang et al.
have used this dataset as a benchmark to evaluate a number of docking scoring functions  <cit> , we used it to compare the performance of motifscore with that of a number of other scoring functions.
atom type assignment
based on the work of others  <cit>  and the chemical properties of simple molecules, we created an atom type classification scheme to describe protein-ligand interactions.
we defined a total of  <dig> atom types , of which  <dig> were for protein atoms and  <dig> for ligand atoms, with many of them shared by both.
the atom types were identified by a 3-code name, or a 3-letter name for those that did not need to be further classified, as they were relatively rarely observed in our dataset of protein-ligand complexes .
some general rules for the 3-code names were as follows.
the 1st code was the name of the element  and the 2nd and 3rd code indicated the surroundings and electrostatic properties of the atom.
the 2nd code could be  <dig>   <dig>  r, or l, which, respectively, correspond to sp <dig> or sp <dig> hybridization or inclusion in an aromatic ring or an aliphatic chain.
the 3rd code could be p, n, a, d, b, e, or c, which, respectively, correspond to polar, non-polar, hydrogen bond acceptor, hydrogen bond donor, an atom that can be both a hydrogen bond donor and acceptor or either a hydrogen bond donor or acceptor, or a charged atom.
the 'either ' code was associated with the atom type nre, which was used primarily for the two nitrogen atoms on the imidazole ring of a histidine, as both of the nitrogens can be either protonated  or deprotonated .
for simplicity, the nitrogen of tryptophan, which is much infrequently seen than histidine especially in the active site, is also assigned the atom type nre.
a metal ions, such as zn, were considered in this work as an extended interaction site of the protein molecule.
b see fig.
s <dig> in additional file  <dig> for a representative chemical structure of the atom type crn, crp, ora, nra, nrd, nre, and sra.
in our scheme, each individual atom is itself an interacting center.
for example, the chemical group po <dig> would have four interacting centers: namely, one pho  and three olc .
others have used pseudoatoms or pseudocenters , such as the geometric center of a chemical group, which may have the advantage of the resulting scoring function being less sensitive to the exact interaction distances used.
on the other hand, chemical group centers do not fully account for the complex interactions, including the different interacting orientations, for example, of the constituent atoms.
the interaction network motifs derived here represent a hybrid of both, since interactions are computed on real atoms but the scoring function is a non-energy-based network motif count, which does not rely on a precise range of interaction distances .
the interaction networks could also be constructed using pseudoatoms instead, but whether this would yield better results remains to be investigated.
for an automated atom type assignment, for protein molecules, we relied on standard pdb files, which have conventional names and thus the implicit bonding information for all protein atoms, and, for ligand molecules, we relied on the het dictionary file on the pdb website  <cit> , which contains the complete naming and bonding information of all the ligand molecules that appear in the pdb files.
determination of atom-atom interaction thresholds
whether or not there is an interaction between two atoms is determined by the energy produced by the interaction, which is affected by many variables, but mostly the distance between the two atoms.
to save computational time, most studies have used distance instead of energy as the criterion to determine whether two atoms interact.
the threshold of distance for interaction could greatly affect the complexity and outcome of a study, but most knowledge-based docking studies seem to suggest that a cutoff in the range of 4- <dig> Å between heavy atoms can achieve optimal performance  <cit> .
however, a single cutoff may not be sufficient for the present study, because: 1) whereas conventional methods use a distance cutoff to avoid the computation of the less significant portion of interaction energies, we used the distance cutoff merely to establish network edges, which, unlike conventional scoring functions, do not involve a distance-dependent energy computation; 2) in the absence of a distance-dependent energy function, clashes between two interacting atoms cannot be prevented.
in order to appropriately define atom-atom interactions , we therefore introduced not only an upper distance threshold for any pair of interacting atoms, but also a lower one.
as described below, the values of these thresholds were determined by examining the distance distributions of all protein-ligand atom pairs  in the  <dig>  protein-ligand complex structures selected from the pdb.
following the concept of radial distribution  <cit> , we analyzed the occurrence distribution of all atom type pairs as a function of separation distance.
for atom type x, in order to determine the distance interval at which the highest density of atom type pair x-y occurs , we divided the occurrence of y by the difference between the cubes of the outer and inner radial distances from atom x .
for normalization, we also divided the occurrence of y atoms found in this particular distance interval by the total occurrence of y atoms in the  <dig>  complex structures in the pdb.
the resulting distributions were quite distinctive from one atom type pair to another.
however, to a large extent, they could all be characterized, albeit approximately, as one of four categories, as summarized in fig.
<dig>
the first category  featured a sharp peak at a small distance, indicating a specific distance or a narrow range of distances for contact  between the two atom types.
more than one-third of the atom type pairs  were classified in this category .
all of these are known to exhibit specific interactions, such as hydrogen bonding, salt bridge, and polar-polar interactions.
examples of this category include hydrogen bonding pairs, such as nlb-o2a and nlb-olc, which represent a hydrogen bonding interaction between a primary or secondary amine  and a double-bonded oxygen  or charged oxygen .
for the second category , the distribution curve increased rapidly over a short span of distance, then stayed fairly flat or decreased slowly over a wide distance range.
this type of distribution is indicative of the absence of a strongly preferred distance for the two atom types to contact each other.
about half of the atom type pairs  fell into this category.
most were for weak interactions , such as c3n-c3n  and crn-crn , or for less well-defined interactions  involving sulfuric or aromatic atom types, such as s3n-crp  and nrd-crn .
the distribution curve of the third category  also began to rise after a short distance, but, unlike the first two categories, rose quite slowly over a large distance.
this indicates that the atom type pairs belonging to this category rarely see each other at short distances.
the few occurrences observed at short distances are likely to be a result of an attraction owing to adjacent atoms that are chemical-bonded to the atom in question.
thirty-eight atom type pairs were classified into this category, which typically exhibited an interaction between a non-polar atom and a polar atom or between two polar atom types with the same polarity, such as c2n-c2p  and nlc-nrd .
the fourth category  showed a noisy curve that did not display an easily identifiable pattern, which may indicate insufficient data, as they usually contained an atom type of rare occurrence, such as halogen in hal-nlb and metal in nlc-met.
there were  <dig> atom type pairs in this category.
a common observation for the first three categories of atom type pairs was a zero occurrence when the distance was short enough, suggesting that a distance threshold could be established and that only above this could attraction overcome repulsion, allowing the pair to be observed.
in addition, the occurrence distribution curves of the first two categories fell after reaching a peak, displaying a specific upper threshold.
the third category did not need an upper threshold, since its atom type pairs preferred not to contact each other.
for the fourth category, we decided to use 4Å and 2Å as the upper and lower threshold values, which are arbitrary choices, but are within the range determined for the other three categories of interaction.
for other rarely observed atoms, such as boron, that are not defined in our simple set of  <dig> atom types , the arbitrarily chosen thresholds of the fourth category along with the associated network motifs can be applied to them to facilitate subsequent computations.
despite the common main features, these occurrence distribution curves were not smooth, making an automated determination of thresholds difficult, so we manually inspected the data and recorded the lower and upper thresholds for each of the  <dig> atom type pairs.
in all, the upper thresholds thus determined ranged from 3Å to  <dig> Å and the lower thresholds from  <dig> Å to 4Å.
.
construction of protein-ligand interaction networks
with the upper and lower thresholds for all the interactions determined, we then constructed protein-ligand interaction networks using the following procedures: 1) carrying out atom type assignment on the atoms of the protein and ligand , 2) checking all the protein-ligand interatomic distances against the distance thresholds of the corresponding atom type pairs, and 3) connecting those that met the thresholds and ignoring those that did not.
as an example, fig.
<dig> shows a 2d representation of the 3d protein-ligand interaction network constructed for the complex of carboxypeptidase a and l-benzylsuccinate .
motif searching and development of motifscore
the next step was to identify motifs from the total of  <dig>  protein-ligand interaction networks constructed as described above.
while a number of algorithms for discovering network motifs have been reported  <cit> , they could not be easily adopted here because the protein-ligand interaction networks constructed in this study are a bipartite network, which is different from the conventional gene-gene or protein-protein interaction networks, in which there is usually only one type of node and for which the existing motif identification algorithms have usually been developed.
moreover, since there were only  <dig> atom types, a number much smaller than the number of interacting atoms, many nodes within the same protein-ligand network were indistinguishable, which is quite different from the situation in other biological networks, where a node usually does not occur more than once in the same network.
in the total of  <dig>  protein-ligand interaction networks constructed, the largest number of protein atoms connected to a ligand atom was  <dig>  thus, a network motif involving just one ligand atom would be those connecting  <dig>  to  <dig> ,  <dig> to  <dig> ..., up to  <dig> to  <dig>  however, mathematically, of these " <dig> interacting with n" motifs, the number of those connecting to  <dig> or  <dig> protein atoms would be the largest .
in the actual data, the  <dig> interacting with  <dig> and  <dig> interacting with  <dig> motifs indeed dominated the  <dig> interacting with n motifs , and, because the number of  <dig> interacting with  <dig> and  <dig> interacting with  <dig> motifs was so huge, the other types of motifs became negligible and a scoring function based on these would be severely biased to the  <dig> interacting with  <dig> and  <dig> interacting with  <dig> motifs.
consequently, we decided to search for other simple motifs, such as those of  <dig> ligand atoms interacting with  <dig> protein atoms, where the scoring function would not be severely biased to a particular type of motif.
there are four distinctive topologies of the  <dig> ligand atoms interacting with  <dig> protein atoms motifs, as shown in fig.
<dig>  in all, we obtained  <dig>  such interaction motifs, the large number arising from the numerous combinations of the  <dig> atom types  that could be assigned to the five constituent atoms.
for various reasons, other simple motifs were not included in motifscore: the  <dig> interacting with  <dig> and  <dig> interacting with  <dig> motifs are likely to occur randomly, as they contain only two connections; there was a comparatively very small number of  <dig>  interacting with  <dig>   <dig> interacting with  <dig>  and  <dig> interacting with  <dig> motifs due to the small number of ligand atoms that can simultaneously interact with protein atoms under our interaction criteria; and, for the  <dig> interacting with  <dig> and  <dig> interacting with  <dig> motifs, no significant difference in the performance of motifscore was seen when they were included .
intuitively, without the guidance of energy computation, a network motif-based scoring function stipulates that the presence of more motifs implicates better interactions .
however, in developing motifscore, we immediately faced two major difficulties: 1) as mentioned, due to the large number of combinations of atom types, even the simple 2- <dig> network topologies  harbored a large number of distinctive motifs which should not all be counted equally, because, for instance, some motifs consist of more frequently occurring atom types than others; 2) if the ligand molecule was pushed into the binding site to make more contacts with the protein, new motifs would be created and some of the old motifs removed, but the rate of increase would far outweigh the rate of loss.
to overcome these difficulties, motifscore was consequently made a composite of two opposing components, gain and penalty, with gain being a sum of the normalized motif counts and penalty being a factor using number of clashes between ligand and protein atoms to avoid motif overcounting, as described below.
in eq.
<dig> we define, for each of the total  <dig>  motifs, a significance grade, sg, which reflects the relative importance of a specific motif in the network.
where omi denotes the occurrence of a specific motif i in the  <dig>  protein-ligand interaction networks constructed from the complex structures in the pdb; motif i is one of  <dig>  motifs, each of which consists of  <dig> atoms,  <dig> on the ligand and  <dig> on the protein, and of 4- <dig> network edges, depending on which of the four motif types  motif i belongs to; ni is the number of motifs having the same motif type as motif i; oaj denotes the occurrence of atom type j  in the  <dig>  networks, where atom type j  is the atom type assigned to one of the  <dig> atoms constituting motif i; and mj is either ml or mp , where ml  is the total number of ligand  atoms.
using normalization by the total occurrence of each motif type and the relative occurrence of each atom type, eq.
<dig> takes into account the fact that motifs consisting of fewer edges or more prevalent atom types are likely to occur more often.
with eq.
<dig>  a motif's sg would not necessarily increase in a larger dataset of protein-ligand complex structures in which motif occurrences are bound to increase.
because each motif's sg reflects the probability of the motif occurring in the interaction network, the probability of it occurring simultaneously with other motifs should be a multiplication of all their sgs.
for easier computation, motif sg was log transformed to convert multiplication into addition.
to score a specific protein-ligand interaction network, i.e.
a specific protein-ligand complex conformation, the gain component  of motifscore is the sum of the motif sg  for all the motifs found in this specific network, or the particular protein-ligand binding conformation from which the network resulted.
where sgk is the motif significance grade  for the kth motif found in the specific network constructed from a specific docking solution, and n is the total number of motifs present in this network.
as mentioned earlier, we need a penalty to counteract the excessive motif gains of a growing interaction network that may occur in docking when one places the ligand molecule closer and closer to its protein receptor.
in eq.
<dig>  nc denotes the total number of clashes between ligand and protein atoms, where one clash is counted when the distance between two atoms, one each from the protein and ligand, is less than the lower threshold of their atom type pair  and weight w is a parameter to be optimized against a subset of lpdb, which served as a training set, as described above.
the value of  <dig> for w appeared to yield the best result, though the results were not very sensitive to the value of this parameter .
a second parameter was also introduced to adjust the lower distance thresholds of the interacting atom type pairs for best fit of the training data.
as described earlier, the lower thresholds were determined based on the observed distribution of interatomic distances.
however, at a distance barely exceeding the lower threshold, the probability of observing the atom type pairs in the native structures is rather low, whereas, when evaluating the conformations produced in the docking process, for any distance that was above the lower threshold, even just barely, a connection/interaction would be made and it would contribute to motif gain.
in order to reflect this reality better, we modified the lower distance thresholds by a factor and found that better results could be obtained when the lower distance thresholds for computing the gain component were raised by  <dig> percent compared to those used to compute penalties , though, like the penalty weight w , the results were not very sensitive to the value of the lower threshold factor used .
finally, after experimenting with several different composite formulae, we settled on one with a ratio for motifscore , which yielded the best results.
RESULTS
scoring docking solutions
as in many related studies, we evaluated the performance of motifscore using the criteria of distance rmsd  between the experimentally observed ligand positions and those of the highest-scored docking solution to consider whether the scoring was a success or a failure.
we used two docking datasets, one by brooks et al.
<cit>  for parameter optimization  and the other by wang et al.
<cit>  for performance comparison with  <dig> other scoring functions .
our success rate using the criterion of rmsd ≦ <dig> Å was 84% for both the parameter optimization set  and wang's benchmark test set .
as shown in table  <dig>  compared to other scoring functions, motifscore performed admirably, being second only to drugscorecsd  <cit> .
drugscorecsd is a significantly improved version of drugscorepdb  <cit> .
drugscorepdb was developed based on the pdb and drugscorecsd was developed based on the cambridge structural database   <cit> , which is a higher-resolution data source for contact distances between interacting atoms, a difference thought to underlie the improvement  <cit> .
whether motifscore can be similarly improved using the csd instead of the pdb remains to be determined, but the effect may not be as significant as in drugscorecsd, as it is likely that the network motif-based scoring  cannot be fine-tuned as much as a conventional energy-based scoring function, such as that used by drugscorecsd.
nevertheless, it is encouraging that, despite its present crude form, motifscore performed surprisingly well and therefore can serve as a non-conventional alternative to existing scoring functions, being useful especially for coarse-grained docking computations.
motifscore
a the data were taken from wang et al.
<cit>  and velec et al.
<cit> , except for the motifscore results, which were computed in this work using the same dataset.
bscoring functions are ranked by their success rates at rmsd ≦  <dig>  Å
cnot provided by velec et al.
<cit>
the 100-complex dataset of wang et al.
<cit>  had also been divided into subsets of hydrophilic, mixed, and hydrophobic complexes to evaluate potential bias of scoring functions on different types of molecular interactions.
as summarized in table  <dig>  motifscore achieved a very high success rate of 91% on both the hydrophilic and mixed subsets, but, like many other scoring functions, including drugscorepdb , its performance was relatively poor for the hydrophobic subset.
the less favorable result for the hydrophobic subset is thought to be due to several difficult cases in this subset in which the binding site is shallow and hard to score accurately  <cit> , compounded by the fact that the number  of complexes in this subset is much smaller than in the other two subsets, which means that the success rate would drop by as much as 4% if the number of correctly predicted complexes was reduced by only one.
the data were taken from wang et al.
<cit> , except for motifscore.
the success rates were determined using rmsd ≦  <dig>  Å.
some significant motifs
to elucidate which of the hundreds of thousands of motifs make the greatest contribution to motifscore, we ranked them by their motif sg.
the top  <dig> motifs are listed in table  <dig>  all contain aromatic and/or polar  atom types on the protein side.
this agrees well with the observation that the frequently observed catalytic residues of proteins are generally the residues with aromatic, polar, or charged side-chains  <cit> .
in addition, although cysteine occurs less often than other amino acids in proteins, it is relatively enriched in catalytic residues  <cit> .
the sulfur atom of cysteines, assigned as s3n, was also commonly observed in our top-ranked motifs.
similarly, there is a relatively small number of halogens in protein binding ligands, but, when present, these atoms usually play a key role in protein-ligand interactions.
the normalized motif count, or sg, reflects the significance of these halogens: atom type hal-containing motifs appear several times in the top  <dig> motifs even though their occurrences are small.
the occurrences were counted in the  <dig>  protein-ligand complex structures used to construct the networks .
see fig.
<dig> for the different motif types.
fig.
4a-d illustrate four of the top  <dig> ranked motifs observed in the  <dig>  protein-ligand complexes.
the four examples all exhibit interactions between aromatic atom types.
as can be seen, these motifs are a composite of multiple pairwise interactions between two ligand interacting sites  and three protein atoms from one , two , or three  amino acid residues.
interestingly, in fig.
4b-d, the aromatic rings of the ligand and of the amino acid side-chains stack against each other, formingπ-πinteractions  <cit> .
these ring-stacking interactions are not easily modeled by conventional scoring functions using separate accounts of pair-wise interactions, but, collectively, they emerge as significant motifs in protein-ligand interaction networks.
binding site-enriched protein triangles
since all of the protein-ligand interaction motifs in our model contain  <dig> protein atoms, we wondered whether the chance of observing three spatially close protein atoms simultaneously was quite different in the ligand-binding site compared to the rest of the protein.
to answer this question, we computed a binding site enrichment factor, fb, using eq.
<dig>  for every distinguishable group of three protein atom types captured in our  <dig> interacting with  <dig> ligand-protein interaction motifs.
for convenience, we called the three protein atoms a protein triangle.
note that, as protein triangles are distinguished by the atom types assigned to their three constituent atoms, they are atom type triangles.
eq.
<dig> was adopted here based on the work of zeeberg et.
al.
<cit> , where ns and ns are the occurrences of a specific atom type triangle found, in, respectively, the interaction network and the whole protein and n and n are the corresponding total number of atom type triangles.
to count the triangles not only in the binding site, but also in the whole protein, we searched for those with all three sides longer than  <dig> Å and two less than  <dig> Å and the other less than  <dig> Å, which would account for almost all  of the triangles present in the protein-ligand interaction motifs identified.
theoretically, with  <dig> protein atom types and taking repetition into account, we should have a total of  <dig> different atom type triangles.
however, as metal ions rarely occur in our dataset of protein structures and almost all of those that do are involved in ligand binding, the fbs for metal-containing triangles were extremely large.
they were therefore excluded in the discussion below, leaving  <dig> non-metal triangles to be ranked and sorted by their fb .
the binding site enrichment factor fb represents a protein triangle's propensity for occurring in the protein-ligand interaction networks constructed, or, in other words, the ligand-binding sites.
as shown in fig.
<dig>  such propensity is far from uniform, and, in fact, the distribution resembles that observed for a wide spectrum of biological properties: the distribution of fb tended to follow a power law distribution  <cit> , y~x-γ, suggesting that only a handful of triangles were highly enriched in ligand-binding sites, while most showed little or no propensity .
some low frequency triangles appear to be concentrated at ligand-binding sites.
for example, the nre-nre-nre triangles, which are triangles that connect  <dig> or  <dig> different histidines , were not commonly observed, but, when they were, they tended to be at ligand-binding sites .
the uneven distribution of fb , which was normalized to account for the wide range of occurrences of these triangles , suggests that fb does not necessarily correlate with binding site occurrence , although triangles with a high fb tended to occur less frequently in other parts of the protein .
many of the highly enriched triangles were constituents of the top-ranked motifs  that were major contributors to the motifscore.
consequently, as can be seen from fig.
<dig>  triangles with higher fb values also tended to form motifs with a higher sg ; the correlation between the two was significant, with the spearman correlation coefficient being  <dig> .
interestingly, table  <dig> shows that some of the binding-site enriched protein triangles consisted of  <dig> charged or polar atom types with the same positive  or negative  charge.
for example, the most enriched triangle consisted of  <dig> nlcs, which are nitrogen atoms with a positive charge on lysine or arginine residues.
the reason why these positively charged amino acids, which are spatially close, as they form the triangle, are enriched in the active site is probably because they all interact with the ligand molecule at the same time.
fb
*coverage: the percentage of protein-ligand complexes  in which a specific triangle was present in the ligand-binding site or the whole protein.
as shown in table  <dig>  the top  <dig> binding site-enriched triangles were generally composed of atom types  that exist on the side-chains of the amino acids histidine, lysine, arginine, asparagine, aspartic acid, glutamic acid, glutamine, serine, tyrosine, threonine, and cysteine.
a survey by barlett et.
al.
<cit>  found that, of the  <dig> amino acids, only these  <dig> polar and charged residues are directly involved in catalytic reactions.
gutteridege and thornton  <cit>  further proposed that histidine is the most commonly observed and most important residue in protein enzymatic reactions, followed by, in descending order of observance, cysteine, the charged residues glutamate, aspartate, arginine, and lysine, and, finally, the polar residues serine, threonine, tyrosine, glutamine, and asparagine.
this is fairly consistent with our statistics on the ligand-binding site enrichment: of the  <dig> protein triangles that were most enriched in ligand-binding sites, crp  was the most commonly observed atom type, while nlc , nre , and olc  were also frequently observed, and nlb  and o3b  on polar residues were observed less often.
unlike in table  <dig>  s3n  appeared only once in table  <dig>  this can be explained by the fact that whereas every protein triangle in table  <dig> is unique, the same triangle can appear multiple times in table  <dig> together with different interacting ligand atoms and/or in different motif topologies .
many of the binding site-enriched triangles are composites of multiple residues.
for example, many of the nre-nre-olc triangles were composed of two histidines and a carboxylate residue, corresponding to the 2-his-1-carboxylate facial triad, a characteristic motif of a metalloenzyme superfamily  <cit> .
similarly, many of the nre-nre-nre triangles were composed of three histidines, known as the histidine triad of the nucleotide-binding histidine triad superfamily  <cit> .
interestingly, while only  <dig> and  <dig> protein structures were annotated by pdb  <cit>  as having, respectively, a histidine-triad or a 2-his-1-carboxylate facial triad in the  <dig>  protein-ligand complex structures we analyzed, roughly 2% and 5% of these proteins contain nre-nre-nre  and nre-nre-olc  atom type triangles in their ligand-binding sites, respectively.
although these proteins may not be members of the histidine triad superfamily or the 2-his-1-carboxylate facial triad superfamily per se, the much higher prevalence than previously recognized of these motifs may guide future investigations to identify conserved catalytic mechanisms in diverse enzyme families.
issues and prospects of motifscore
although motifscore has been derived using static crystallographic structures of protein-ligand complexes, the use of non-precise interaction distances and counts of number of motifs instead of interaction energies likely renders it more tolerant to subtle conformational changes in protein or ligand than are conventional energy-based scoring functions, which are known to be sensitive to steric clashes especially in unbound dockings  <cit> .
the results on both the lpdb and wang's dataset attest the ability of motifscore to account for at least ligand flexibilities, since the ligands, though not the proteins, in the decoys of both sets include flexible conformations generated by molecular dynamics simulation and genetic algorithm, respectively.
nevertheless, the performance of motifscore in dockings that sample different protein conformations  requires further examinations.
motifscore can be regarded as a kind of knowledge-based scoring function since it was also derived by extracting information from a statistical analysis of known protein-ligand complex structures.
however, motifscore is different from conventional knowledge-based scoring functions such as pmf and drugscore  <cit>  in at least two aspects: 1) it is non-energy based, and 2) using the interaction network motifs, it can score directly on the 3d interaction patterns of molecular recognition conserved in protein ligand complexes.
one disadvantage of motifscore, however, is that a non-conventional search scheme may need to be developed to take advantages of its unique features, as discussed next.
as motifscore is non-energy-based, it was of interest to examine the landscape of its functional values.
we found that motifscore did not correlate well with experimentally determined binding affinities, the spearmen correlation coefficient between the two being only  <dig> , only better than that  of autodock among the various scoring functions evaluated by wang et al.
<cit> .
nevertheless, fig.
<dig> shows that, for a typical success case, motifscore could easily distinguish reasonably good docking solutions  from bad ones.
intriguingly, there appears to be a very narrow funnel leading to the native state formed by very good docking solutions.
this is quite distinct from that of a conventional energy-based scoring function where the funnel leading to the energy minimum, which should be reasonably close to the native state, is usually much smoother  <cit> .
this implies that, whereas a search algorithm, such as a genetic algorithm, may work efficiently with a conventional, energy-based scoring function to find good docking solutions, its direct adoption for use with motifscore would not be ideal.
on the other hand, the protein-ligand interaction motifs derived in this work, in their three dimensional arrangements of interacting atoms, capture the spatial arrangements of reasonably good docking solutions, so it is quite possible to develop a scheme to look up the table of interaction network motifs  and directly home in on a reasonable docking solution based on the structure of a few top-ranked motifs.
this could eliminate the time-consuming searching computation needed in conventional docking methods, or at least provide a good starting solution for further refinement.
in addition, although motifscore is currently limited to scoring the interactions between protein and small molecules, the same network-based approach should be extendable to protein-protein and protein-dna/rna interactions, which form conserved spatial and chemical binding patterns .
work along these lines is currently undergoing in our laboratory.
CONCLUSIONS
motifscore is a novel interaction-motif-based scoring function for protein-ligand docking.
despite the absence of mathematical models to mimic the force field of molecular interactions, motifscore performed well in distinguishing between good and bad docking solutions in a benchmark test set.
furthermore, owing to the network approach, motifscore is intrinsically more able than conventional docking scoring functions to capture interactions involving more than two interacting sites, the π-π stacking of two aromatic rings being a prime example.
the ligand-binding site-enriched interaction motifs identified are in accord with existing knowledge on protein-ligand binding and may prove useful for binding site predictions.
finally, the three-dimensional protein-ligand interacting motifs could provide very good templates for placing ligand molecules in fast, though coarse, protein docking computations.
availability and requirements
the source code of motifscore is available for download at: https://sourceforge.net/projects/msdock/.
the package contains a set of perl scripts for computing the scoring function and for creating a perl database file that stores the names of ligands and their atom types.
it also offers some demonstration examples of how to obtain the final docking scores.
these scripts work on a unix or linux platform and their download is free for academic users.
authors' contributions
zrx carried out the research and drafted the manuscript.
mjh conceived the project, guided the research, and revised the manuscript.
both authors read and approved the final manuscript.
supplementary material
additional file 1
supplementary data for details of motif related parameters.
this file contains  <dig> tables and  <dig> figures.
table s <dig> lists the lower and upper distance thresholds for each of the  <dig> atom type pairs.
table s <dig> lists the binding site enrichment factor, fb, of each of the  <dig> non-metal protein atom type triangles.
figs.
s <dig> and s <dig> display the optimization results of two parameters, penalty weight w and the raising percentage of lower thresholds, on the training dataset lpdb.
fig.
s <dig> shows a schematic presentation for some of the atom types defined in table  <dig>
click here for file
acknowledgements
we would like to thank dr. shaomeng wang  and dr. charles l. brooks iii  for making their docking datasets available, and dr. jung-hsin lin  and dr. yen-jen oyang  for helpful discussions.
the work was supported by the national science council  .

 drug design and discovery play a pivotal role in driving research in computational chemistry and biology. in computational drug design and discovery, it is often necessary to determine, as a first step, the binding of a ligand to a target protein. the computational scheme for predicting ligand binding occurrence, affinity, and orientation is commonly referred to as "molecular docking", which has been a topic of intensive research for decades. the development of a molecular docking tool usually starts with an efficient search algorithm, which places the ligand in the active site of the target protein in numerous different positions, orientations, and, in flexible docking, conformations. these are then evaluated by a scoring function to distinguish between good and bad docking solutions. the two aspects of searching and scoring can be, and usually have been, developed and evaluated separately, although one clearly affects the other and a balance is often sought to meet specific study aims. scoring functions for molecular docking are traditionally either physics-based or knowledge-based and differ mainly in the derivation of the mathematical models used to compute the energies of molecular interactions. physics-based scoring functions employ a set of molecular interaction terms to compute binding energies. for example, the scoring function used in g-score or autodock is based on the molecular mechanics force field used in tripos or amber, respectively, and f-score and chemscore derive the coefficients of their energy terms using regression analysis of experimentally determined binding energies. in contrast, knowledge-based scoring functions, such as pmf and drugscore, rely on statistical observations of preferred protein-ligand contacts, from which binding energies are calculated. a common thread in these scoring functions is that, with few exceptions, they all operate on the assumption that the total interaction will be faithfully represented by an additive summation over a series of pairwise interactions between interaction centers, which are either real or imaginary. despite fundamental fallacies in this assumption, it has been adopted for many years, as otherwise the problem would be intractable or too computationally expensive to handle practically. in this study, we investigated a new idea for developing a scoring function for molecular docking. the novelty of the idea lies in the use of a network approach to identify frequently occurring patterns of protein-ligand interactions not in the form of pairwise interactions, but in the form of network motifs, where "network" refers to a collection of pairwise interactions. the motifs, which are assemblies of multiple pairwise interactions between proteins and ligands that are frequently observed in a database of known protein-ligand complex structures, thus represent energetically favorable ways of positioning a ligand molecule in the active site of a protein. such a motif-based approach offers a new line of exploration for research on molecular docking. as a first step in this line of research, we present the construction of a motif-based scoring function, called motifscore, and the evaluation of its ability to distinguish between good and bad docking solutions. our results using a benchmark dataset showed that motifscore performed well against a number of scoring functions used in existing docking programs. methods in constructing the protein-ligand interaction network, the three-dimensional coordinates of a set of protein-ligand complexes were each transformed into an atom-atom interaction network. there are two types of nodes in the network, namely atoms from the protein and the ligand. two nodes are connected by an edge if their interaction is deemed significant by a distance threshold determined from a statistical analysis of a dataset of protein-ligand complexes. interactions between atoms within the protein or within the ligand are not considered. these networks were then broken down into many interaction motifs representing simple units of the network with specific protein-ligand interaction patterns. datasets to establish distance thresholds for constructing the protein-ligand interaction network, we performed a large-scale statistical analysis on a diverse set of known structures of protein-ligand complexes. the available protein-ligand complex structures as of january in the protein data bank were screened. to simplify our task, only those determined by x-ray crystallography and with only one ligand molecule were selected. furthermore, those containing dna/rna molecules and those in which the ligand binds to its target protein covalently were discarded. complexes containing a heme group were also excluded because one may consider the heme group as a part of the protein, rather than a ligand. solvent molecules and ions, such as chloride or ammonium, that are often included to facilitate crystallographic studies were not considered as ligands, whereas metal ions such as zinc or calcium ion that are coordinated by metal-binding amino acids were considered as an extension of the protein molecule, and their interactions with ligand were considered part of the protein-ligand interaction network constructed. in all, we obtained a total of structures of complexes that could be used in the statistical analysis. this set of pdb protein-ligand complex structures is not entirely non-redundant since a single protein or its mutants may be co-crystallized with many similar compounds, and, conversely, the same compound may be co-crystallized with multiple homologous proteins. however, unless both the protein and the ligand are completely identical, part of the resulting interaction networks can still be distinct, and as defining non-redundant networks is not a straightforward task, we chose to include as many network connectivities as possible while conducting normalizations to minimize the potential bias of using a redundant dataset. furthermore, analysis on a reduced set by excluding those that share >25% protein sequence identity and identical ligand showed that the resulting network motifs and their grades had a pearson's correlation coefficient of with those of the set of structures, suggesting that including homologous proteins and identical ligand compounds would not significantly affect the outcome of our analysis. besides the set of structures, we used a subset of the lpdb, which contains decoys generated from molecular dynamics simulations, to optimize the parameters of our scoring function, motifscore. there are more than complexes, each with a set of various numbers of decoys, in the lpdb, but only of these were used in this work because we found that the ligand names, atom names, and the ordering of atoms of many of the decoys in lpdb were not the same as in the original data files in the pdb, making it difficult for us to correctly assign their designated atom types. for testing and evaluating motifscore, the data set of wang et al. was used. the wang dataset consists of a set of protein-ligand complexes, each of which comes with docked conformations generated using autodock. since wang et al. have used this dataset as a benchmark to evaluate a number of docking scoring functions, we used it to compare the performance of motifscore with that of a number of other scoring functions. atom type assignment based on the work of others and the chemical properties of simple molecules, we created an atom type classification scheme to describe protein-ligand interactions. we defined a total of atom types, of which were for protein atoms and for ligand atoms, with many of them shared by both. the atom types were identified by a 3-code name, or a 3-letter name for those that did not need to be further classified, as they were relatively rarely observed in our dataset of protein-ligand complexes. some general rules for the 3-code names were as follows. the 1st code was the name of the element and the 2nd and 3rd code indicated the surroundings and electrostatic properties of the atom. the 2nd code could be r, or l, which, respectively, correspond to sp or sp hybridization or inclusion in an aromatic ring or an aliphatic chain. the 3rd code could be p, n, a, d, b, e, or c, which, respectively, correspond to polar, non-polar, hydrogen bond acceptor, hydrogen bond donor, an atom that can be both a hydrogen bond donor and acceptor or either a hydrogen bond donor or acceptor, or a charged atom. the'either' code was associated with the atom type nre, which was used primarily for the two nitrogen atoms on the imidazole ring of a histidine, as both of the nitrogens can be either protonated or deprotonated. for simplicity, the nitrogen of tryptophan, which is much infrequently seen than histidine especially in the active site, is also assigned the atom type nre. a metal ions, such as zn, were considered in this work as an extended interaction site of the protein molecule. b see fig. s in additional file for a representative chemical structure of the atom type crn, crp, ora, nra, nrd, nre, and sra. in our scheme, each individual atom is itself an interacting center. for example, the chemical group po would have four interacting centers: namely, one pho and three olc. others have used pseudoatoms or pseudocenters, such as the geometric center of a chemical group, which may have the advantage of the resulting scoring function being less sensitive to the exact interaction distances used. on the other hand, chemical group centers do not fully account for the complex interactions, including the different interacting orientations, for example, of the constituent atoms. the interaction network motifs derived here represent a hybrid of both, since interactions are computed on real atoms but the scoring function is a non-energy-based network motif count, which does not rely on a precise range of interaction distances. the interaction networks could also be constructed using pseudoatoms instead, but whether this would yield better results remains to be investigated. for an automated atom type assignment, for protein molecules, we relied on standard pdb files, which have conventional names and thus the implicit bonding information for all protein atoms, and, for ligand molecules, we relied on the het dictionary file on the pdb website, which contains the complete naming and bonding information of all the ligand molecules that appear in the pdb files. determination of atom-atom interaction thresholds whether or not there is an interaction between two atoms is determined by the energy produced by the interaction, which is affected by many variables, but mostly the distance between the two atoms. to save computational time, most studies have used distance instead of energy as the criterion to determine whether two atoms interact. the threshold of distance for interaction could greatly affect the complexity and outcome of a study, but most knowledge-based docking studies seem to suggest that a cutoff in the range of 4- Å between heavy atoms can achieve optimal performance. however, a single cutoff may not be sufficient for the present study, because: 1) whereas conventional methods use a distance cutoff to avoid the computation of the less significant portion of interaction energies, we used the distance cutoff merely to establish network edges, which, unlike conventional scoring functions, do not involve a distance-dependent energy computation; 2) in the absence of a distance-dependent energy function, clashes between two interacting atoms cannot be prevented. in order to appropriately define atom-atom interactions, we therefore introduced not only an upper distance threshold for any pair of interacting atoms, but also a lower one. as described below, the values of these thresholds were determined by examining the distance distributions of all protein-ligand atom pairs in the protein-ligand complex structures selected from the pdb. following the concept of radial distribution, we analyzed the occurrence distribution of all atom type pairs as a function of separation distance. for atom type x, in order to determine the distance interval at which the highest density of atom type pair x-y occurs, we divided the occurrence of y by the difference between the cubes of the outer and inner radial distances from atom x. for normalization, we also divided the occurrence of y atoms found in this particular distance interval by the total occurrence of y atoms in the complex structures in the pdb. the resulting distributions were quite distinctive from one atom type pair to another. however, to a large extent, they could all be characterized, albeit approximately, as one of four categories, as summarized in fig. the first category featured a sharp peak at a small distance, indicating a specific distance or a narrow range of distances for contact between the two atom types. more than one-third of the atom type pairs were classified in this category. all of these are known to exhibit specific interactions, such as hydrogen bonding, salt bridge, and polar-polar interactions. examples of this category include hydrogen bonding pairs, such as nlb-o2a and nlb-olc, which represent a hydrogen bonding interaction between a primary or secondary amine and a double-bonded oxygen or charged oxygen. for the second category, the distribution curve increased rapidly over a short span of distance, then stayed fairly flat or decreased slowly over a wide distance range. this type of distribution is indicative of the absence of a strongly preferred distance for the two atom types to contact each other. about half of the atom type pairs fell into this category. most were for weak interactions, such as c3n-c3n and crn-crn, or for less well-defined interactions involving sulfuric or aromatic atom types, such as s3n-crp and nrd-crn. the distribution curve of the third category also began to rise after a short distance, but, unlike the first two categories, rose quite slowly over a large distance. this indicates that the atom type pairs belonging to this category rarely see each other at short distances. the few occurrences observed at short distances are likely to be a result of an attraction owing to adjacent atoms that are chemical-bonded to the atom in question. thirty-eight atom type pairs were classified into this category, which typically exhibited an interaction between a non-polar atom and a polar atom or between two polar atom types with the same polarity, such as c2n-c2p and nlc-nrd. the fourth category showed a noisy curve that did not display an easily identifiable pattern, which may indicate insufficient data, as they usually contained an atom type of rare occurrence, such as halogen in hal-nlb and metal in nlc-met. there were atom type pairs in this category. a common observation for the first three categories of atom type pairs was a zero occurrence when the distance was short enough, suggesting that a distance threshold could be established and that only above this could attraction overcome repulsion, allowing the pair to be observed. in addition, the occurrence distribution curves of the first two categories fell after reaching a peak, displaying a specific upper threshold. the third category did not need an upper threshold, since its atom type pairs preferred not to contact each other. for the fourth category, we decided to use 4Å and 2Å as the upper and lower threshold values, which are arbitrary choices, but are within the range determined for the other three categories of interaction. for other rarely observed atoms, such as boron, that are not defined in our simple set of atom types, the arbitrarily chosen thresholds of the fourth category along with the associated network motifs can be applied to them to facilitate subsequent computations. despite the common main features, these occurrence distribution curves were not smooth, making an automated determination of thresholds difficult, so we manually inspected the data and recorded the lower and upper thresholds for each of the atom type pairs. in all, the upper thresholds thus determined ranged from 3Å to Å and the lower thresholds from Å to 4Å.. construction of protein-ligand interaction networks with the upper and lower thresholds for all the interactions determined, we then constructed protein-ligand interaction networks using the following procedures: 1) carrying out atom type assignment on the atoms of the protein and ligand, 2) checking all the protein-ligand interatomic distances against the distance thresholds of the corresponding atom type pairs, and 3) connecting those that met the thresholds and ignoring those that did not. as an example, fig. shows a 2d representation of the 3d protein-ligand interaction network constructed for the complex of carboxypeptidase a and l-benzylsuccinate. motif searching and development of motifscore the next step was to identify motifs from the total of protein-ligand interaction networks constructed as described above. while a number of algorithms for discovering network motifs have been reported, they could not be easily adopted here because the protein-ligand interaction networks constructed in this study are a bipartite network, which is different from the conventional gene-gene or protein-protein interaction networks, in which there is usually only one type of node and for which the existing motif identification algorithms have usually been developed. moreover, since there were only atom types, a number much smaller than the number of interacting atoms, many nodes within the same protein-ligand network were indistinguishable, which is quite different from the situation in other biological networks, where a node usually does not occur more than once in the same network. in the total of protein-ligand interaction networks constructed, the largest number of protein atoms connected to a ligand atom was thus, a network motif involving just one ligand atom would be those connecting to, to..., up to to however, mathematically, of these " interacting with n" motifs, the number of those connecting to or protein atoms would be the largest. in the actual data, the interacting with and interacting with motifs indeed dominated the interacting with n motifs, and, because the number of interacting with and interacting with motifs was so huge, the other types of motifs became negligible and a scoring function based on these would be severely biased to the interacting with and interacting with motifs. consequently, we decided to search for other simple motifs, such as those of ligand atoms interacting with protein atoms, where the scoring function would not be severely biased to a particular type of motif. there are four distinctive topologies of the ligand atoms interacting with protein atoms motifs, as shown in fig. in all, we obtained such interaction motifs, the large number arising from the numerous combinations of the atom types that could be assigned to the five constituent atoms. for various reasons, other simple motifs were not included in motifscore: the interacting with and interacting with motifs are likely to occur randomly, as they contain only two connections; there was a comparatively very small number of interacting with interacting with and interacting with motifs due to the small number of ligand atoms that can simultaneously interact with protein atoms under our interaction criteria; and, for the interacting with and interacting with motifs, no significant difference in the performance of motifscore was seen when they were included. intuitively, without the guidance of energy computation, a network motif-based scoring function stipulates that the presence of more motifs implicates better interactions. however, in developing motifscore, we immediately faced two major difficulties: 1) as mentioned, due to the large number of combinations of atom types, even the simple 2- network topologies harbored a large number of distinctive motifs which should not all be counted equally, because, for instance, some motifs consist of more frequently occurring atom types than others; 2) if the ligand molecule was pushed into the binding site to make more contacts with the protein, new motifs would be created and some of the old motifs removed, but the rate of increase would far outweigh the rate of loss. to overcome these difficulties, motifscore was consequently made a composite of two opposing components, gain and penalty, with gain being a sum of the normalized motif counts and penalty being a factor using number of clashes between ligand and protein atoms to avoid motif overcounting, as described below. in eq. we define, for each of the total motifs, a significance grade, sg, which reflects the relative importance of a specific motif in the network. where omi denotes the occurrence of a specific motif i in the protein-ligand interaction networks constructed from the complex structures in the pdb; motif i is one of motifs, each of which consists of atoms, on the ligand and on the protein, and of 4- network edges, depending on which of the four motif types motif i belongs to; ni is the number of motifs having the same motif type as motif i; oaj denotes the occurrence of atom type j in the networks, where atom type j is the atom type assigned to one of the atoms constituting motif i; and mj is either ml or mp, where ml is the total number of ligand atoms. using normalization by the total occurrence of each motif type and the relative occurrence of each atom type, eq. takes into account the fact that motifs consisting of fewer edges or more prevalent atom types are likely to occur more often. with eq. a motif's sg would not necessarily increase in a larger dataset of protein-ligand complex structures in which motif occurrences are bound to increase. because each motif's sg reflects the probability of the motif occurring in the interaction network, the probability of it occurring simultaneously with other motifs should be a multiplication of all their sgs. for easier computation, motif sg was log transformed to convert multiplication into addition. to score a specific protein-ligand interaction network, i.e. a specific protein-ligand complex conformation, the gain component of motifscore is the sum of the motif sg for all the motifs found in this specific network, or the particular protein-ligand binding conformation from which the network resulted. where sgk is the motif significance grade for the kth motif found in the specific network constructed from a specific docking solution, and n is the total number of motifs present in this network. as mentioned earlier, we need a penalty to counteract the excessive motif gains of a growing interaction network that may occur in docking when one places the ligand molecule closer and closer to its protein receptor. in eq. nc denotes the total number of clashes between ligand and protein atoms, where one clash is counted when the distance between two atoms, one each from the protein and ligand, is less than the lower threshold of their atom type pair and weight w is a parameter to be optimized against a subset of lpdb, which served as a training set, as described above. the value of for w appeared to yield the best result, though the results were not very sensitive to the value of this parameter. a second parameter was also introduced to adjust the lower distance thresholds of the interacting atom type pairs for best fit of the training data. as described earlier, the lower thresholds were determined based on the observed distribution of interatomic distances. however, at a distance barely exceeding the lower threshold, the probability of observing the atom type pairs in the native structures is rather low, whereas, when evaluating the conformations produced in the docking process, for any distance that was above the lower threshold, even just barely, a connection/interaction would be made and it would contribute to motif gain. in order to reflect this reality better, we modified the lower distance thresholds by a factor and found that better results could be obtained when the lower distance thresholds for computing the gain component were raised by percent compared to those used to compute penalties, though, like the penalty weight w, the results were not very sensitive to the value of the lower threshold factor used. finally, after experimenting with several different composite formulae, we settled on one with a ratio for motifscore, which yielded the best results.  scoring docking solutions as in many related studies, we evaluated the performance of motifscore using the criteria of distance rmsd between the experimentally observed ligand positions and those of the highest-scored docking solution to consider whether the scoring was a success or a failure. we used two docking datasets, one by brooks et al. for parameter optimization and the other by wang et al. for performance comparison with other scoring functions. our success rate using the criterion of rmsd ≦ Å was 84% for both the parameter optimization set and wang's benchmark test set. as shown in table compared to other scoring functions, motifscore performed admirably, being second only to drugscorecsd. drugscorecsd is a significantly improved version of drugscorepdb. drugscorepdb was developed based on the pdb and drugscorecsd was developed based on the cambridge structural database, which is a higher-resolution data source for contact distances between interacting atoms, a difference thought to underlie the improvement. whether motifscore can be similarly improved using the csd instead of the pdb remains to be determined, but the effect may not be as significant as in drugscorecsd, as it is likely that the network motif-based scoring cannot be fine-tuned as much as a conventional energy-based scoring function, such as that used by drugscorecsd. nevertheless, it is encouraging that, despite its present crude form, motifscore performed surprisingly well and therefore can serve as a non-conventional alternative to existing scoring functions, being useful especially for coarse-grained docking computations. motifscore a the data were taken from wang et al. and velec et al., except for the motifscore results, which were computed in this work using the same dataset. bscoring functions are ranked by their success rates at rmsd ≦ Å cnot provided by velec et al. the 100-complex dataset of wang et al. had also been divided into subsets of hydrophilic, mixed, and hydrophobic complexes to evaluate potential bias of scoring functions on different types of molecular interactions. as summarized in table motifscore achieved a very high success rate of 91% on both the hydrophilic and mixed subsets, but, like many other scoring functions, including drugscorepdb, its performance was relatively poor for the hydrophobic subset. the less favorable result for the hydrophobic subset is thought to be due to several difficult cases in this subset in which the binding site is shallow and hard to score accurately, compounded by the fact that the number of complexes in this subset is much smaller than in the other two subsets, which means that the success rate would drop by as much as 4% if the number of correctly predicted complexes was reduced by only one. the data were taken from wang et al., except for motifscore. the success rates were determined using rmsd ≦ Å. some significant motifs to elucidate which of the hundreds of thousands of motifs make the greatest contribution to motifscore, we ranked them by their motif sg. the top motifs are listed in table all contain aromatic and/or polar atom types on the protein side. this agrees well with the observation that the frequently observed catalytic residues of proteins are generally the residues with aromatic, polar, or charged side-chains. in addition, although cysteine occurs less often than other amino acids in proteins, it is relatively enriched in catalytic residues. the sulfur atom of cysteines, assigned as s3n, was also commonly observed in our top-ranked motifs. similarly, there is a relatively small number of halogens in protein binding ligands, but, when present, these atoms usually play a key role in protein-ligand interactions. the normalized motif count, or sg, reflects the significance of these halogens: atom type hal-containing motifs appear several times in the top motifs even though their occurrences are small. the occurrences were counted in the protein-ligand complex structures used to construct the networks. see fig. for the different motif types. fig. 4a-d illustrate four of the top ranked motifs observed in the protein-ligand complexes. the four examples all exhibit interactions between aromatic atom types. as can be seen, these motifs are a composite of multiple pairwise interactions between two ligand interacting sites and three protein atoms from one, two, or three amino acid residues. interestingly, in fig. 4b-d, the aromatic rings of the ligand and of the amino acid side-chains stack against each other, formingπ-πinteractions. these ring-stacking interactions are not easily modeled by conventional scoring functions using separate accounts of pair-wise interactions, but, collectively, they emerge as significant motifs in protein-ligand interaction networks. binding site-enriched protein triangles since all of the protein-ligand interaction motifs in our model contain protein atoms, we wondered whether the chance of observing three spatially close protein atoms simultaneously was quite different in the ligand-binding site compared to the rest of the protein. to answer this question, we computed a binding site enrichment factor, fb, using eq. for every distinguishable group of three protein atom types captured in our interacting with ligand-protein interaction motifs. for convenience, we called the three protein atoms a protein triangle. note that, as protein triangles are distinguished by the atom types assigned to their three constituent atoms, they are atom type triangles. eq. was adopted here based on the work of zeeberg et. al., where ns and ns are the occurrences of a specific atom type triangle found, in, respectively, the interaction network and the whole protein and n and n are the corresponding total number of atom type triangles. to count the triangles not only in the binding site, but also in the whole protein, we searched for those with all three sides longer than Å and two less than Å and the other less than Å, which would account for almost all of the triangles present in the protein-ligand interaction motifs identified. theoretically, with protein atom types and taking repetition into account, we should have a total of different atom type triangles. however, as metal ions rarely occur in our dataset of protein structures and almost all of those that do are involved in ligand binding, the fbs for metal-containing triangles were extremely large. they were therefore excluded in the discussion below, leaving non-metal triangles to be ranked and sorted by their fb. the binding site enrichment factor fb represents a protein triangle's propensity for occurring in the protein-ligand interaction networks constructed, or, in other words, the ligand-binding sites. as shown in fig. such propensity is far from uniform, and, in fact, the distribution resembles that observed for a wide spectrum of biological properties: the distribution of fb tended to follow a power law distribution, y~x-γ, suggesting that only a handful of triangles were highly enriched in ligand-binding sites, while most showed little or no propensity. some low frequency triangles appear to be concentrated at ligand-binding sites. for example, the nre-nre-nre triangles, which are triangles that connect or different histidines, were not commonly observed, but, when they were, they tended to be at ligand-binding sites. the uneven distribution of fb, which was normalized to account for the wide range of occurrences of these triangles, suggests that fb does not necessarily correlate with binding site occurrence, although triangles with a high fb tended to occur less frequently in other parts of the protein. many of the highly enriched triangles were constituents of the top-ranked motifs that were major contributors to the motifscore. consequently, as can be seen from fig. triangles with higher fb values also tended to form motifs with a higher sg; the correlation between the two was significant, with the spearman correlation coefficient being. interestingly, table shows that some of the binding-site enriched protein triangles consisted of charged or polar atom types with the same positive or negative charge. for example, the most enriched triangle consisted of nlcs, which are nitrogen atoms with a positive charge on lysine or arginine residues. the reason why these positively charged amino acids, which are spatially close, as they form the triangle, are enriched in the active site is probably because they all interact with the ligand molecule at the same time. fb *coverage: the percentage of protein-ligand complexes in which a specific triangle was present in the ligand-binding site or the whole protein. as shown in table the top binding site-enriched triangles were generally composed of atom types that exist on the side-chains of the amino acids histidine, lysine, arginine, asparagine, aspartic acid, glutamic acid, glutamine, serine, tyrosine, threonine, and cysteine. a survey by barlett et. al. found that, of the amino acids, only these polar and charged residues are directly involved in catalytic reactions. gutteridege and thornton further proposed that histidine is the most commonly observed and most important residue in protein enzymatic reactions, followed by, in descending order of observance, cysteine, the charged residues glutamate, aspartate, arginine, and lysine, and, finally, the polar residues serine, threonine, tyrosine, glutamine, and asparagine. this is fairly consistent with our statistics on the ligand-binding site enrichment: of the protein triangles that were most enriched in ligand-binding sites, crp was the most commonly observed atom type, while nlc, nre, and olc were also frequently observed, and nlb and o3b on polar residues were observed less often. unlike in table s3n appeared only once in table this can be explained by the fact that whereas every protein triangle in table is unique, the same triangle can appear multiple times in table together with different interacting ligand atoms and/or in different motif topologies. many of the binding site-enriched triangles are composites of multiple residues. for example, many of the nre-nre-olc triangles were composed of two histidines and a carboxylate residue, corresponding to the 2-his-1-carboxylate facial triad, a characteristic motif of a metalloenzyme superfamily. similarly, many of the nre-nre-nre triangles were composed of three histidines, known as the histidine triad of the nucleotide-binding histidine triad superfamily. interestingly, while only and protein structures were annotated by pdb as having, respectively, a histidine-triad or a 2-his-1-carboxylate facial triad in the protein-ligand complex structures we analyzed, roughly 2% and 5% of these proteins contain nre-nre-nre and nre-nre-olc atom type triangles in their ligand-binding sites, respectively. although these proteins may not be members of the histidine triad superfamily or the 2-his-1-carboxylate facial triad superfamily per se, the much higher prevalence than previously recognized of these motifs may guide future investigations to identify conserved catalytic mechanisms in diverse enzyme families. issues and prospects of motifscore although motifscore has been derived using static crystallographic structures of protein-ligand complexes, the use of non-precise interaction distances and counts of number of motifs instead of interaction energies likely renders it more tolerant to subtle conformational changes in protein or ligand than are conventional energy-based scoring functions, which are known to be sensitive to steric clashes especially in unbound dockings. the results on both the lpdb and wang's dataset attest the ability of motifscore to account for at least ligand flexibilities, since the ligands, though not the proteins, in the decoys of both sets include flexible conformations generated by molecular dynamics simulation and genetic algorithm, respectively. nevertheless, the performance of motifscore in dockings that sample different protein conformations requires further examinations. motifscore can be regarded as a kind of knowledge-based scoring function since it was also derived by extracting information from a statistical analysis of known protein-ligand complex structures. however, motifscore is different from conventional knowledge-based scoring functions such as pmf and drugscore in at least two aspects: 1) it is non-energy based, and 2) using the interaction network motifs, it can score directly on the 3d interaction patterns of molecular recognition conserved in protein ligand complexes. one disadvantage of motifscore, however, is that a non-conventional search scheme may need to be developed to take advantages of its unique features, as discussed next. as motifscore is non-energy-based, it was of interest to examine the landscape of its functional values. we found that motifscore did not correlate well with experimentally determined binding affinities, the spearmen correlation coefficient between the two being only, only better than that of autodock among the various scoring functions evaluated by wang et al.. nevertheless, fig. shows that, for a typical success case, motifscore could easily distinguish reasonably good docking solutions from bad ones. intriguingly, there appears to be a very narrow funnel leading to the native state formed by very good docking solutions. this is quite distinct from that of a conventional energy-based scoring function where the funnel leading to the energy minimum, which should be reasonably close to the native state, is usually much smoother. this implies that, whereas a search algorithm, such as a genetic algorithm, may work efficiently with a conventional, energy-based scoring function to find good docking solutions, its direct adoption for use with motifscore would not be ideal. on the other hand, the protein-ligand interaction motifs derived in this work, in their three dimensional arrangements of interacting atoms, capture the spatial arrangements of reasonably good docking solutions, so it is quite possible to develop a scheme to look up the table of interaction network motifs and directly home in on a reasonable docking solution based on the structure of a few top-ranked motifs. this could eliminate the time-consuming searching computation needed in conventional docking methods, or at least provide a good starting solution for further refinement. in addition, although motifscore is currently limited to scoring the interactions between protein and small molecules, the same network-based approach should be extendable to protein-protein and protein-dna/rna interactions, which form conserved spatial and chemical binding patterns. work along these lines is currently undergoing in our laboratory.  motifscore is a novel interaction-motif-based scoring function for protein-ligand docking. despite the absence of mathematical models to mimic the force field of molecular interactions, motifscore performed well in distinguishing between good and bad docking solutions in a benchmark test set. furthermore, owing to the network approach, motifscore is intrinsically more able than conventional docking scoring functions to capture interactions involving more than two interacting sites, the π-π stacking of two aromatic rings being a prime example. the ligand-binding site-enriched interaction motifs identified are in accord with existing knowledge on protein-ligand binding and may prove useful for binding site predictions. finally, the three-dimensional protein-ligand interacting motifs could provide very good templates for placing ligand molecules in fast, though coarse, protein docking computations. availability and requirements the source code of motifscore is available for download at: https://sourceforge.net/projects/msdock/. the package contains a set of perl scripts for computing the scoring function and for creating a perl database file that stores the names of ligands and their atom types. it also offers some demonstration examples of how to obtain the final docking scores. these scripts work on a unix or linux platform and their download is free for academic users. authors' contributions zrx carried out the research and drafted the manuscript. mjh conceived the project, guided the research, and revised the manuscript. both authors read and approved the final manuscript. supplementary material additional file 1 supplementary data for details of motif related parameters. this file contains tables and figures. table s lists the lower and upper distance thresholds for each of the atom type pairs. table s lists the binding site enrichment factor, fb, of each of the non-metal protein atom type triangles. figs. s and s display the optimization results of two parameters, penalty weight w and the raising percentage of lower thresholds, on the training dataset lpdb. fig. s shows a schematic presentation for some of the atom types defined in table click here for file acknowledgements we would like to thank dr. shaomeng wang and dr. charles l. brooks iii for making their docking datasets available, and dr. jung-hsin lin and dr. yen-jen oyang for helpful discussions. the work was supported by the national science council.
it is non-energy-based, and docking is, instead, scored by counting the occurrences of motifs of protein-ligand interaction networks constructed using structures of protein-ligand complexes. in conventional scoring functions, energy terms modeling pairwise interactions are cumulatively summed, and the best docking solution is selected. many of the top motifs, which comprise a multitude of chemical groups that interact simultaneously and make a highly significant contribution to motifscore, capture recurrent interacting patterns beyond pairwise interactions. in this benchmark test, 84% of the highest-scored docking conformations had root-mean-square deviations below Å from the native conformation, which is comparable with the best of several energy-based docking scoring functions. here, we propose to transform protein-ligand interactions into three-dimensional geometric networks, from which recurring network substructures, or network motifs, are selected and used to provide probability-ranked interaction templates with which to score docking solutions.
 a good scoring function is essential for molecular docking computations. in conventional scoring functions, energy terms modeling pairwise interactions are cumulatively summed, and the best docking solution is selected. here, we propose to transform protein-ligand interactions into three-dimensional geometric networks, from which recurring network substructures, or network motifs, are selected and used to provide probability-ranked interaction templates with which to score docking solutions.  a novel scoring function for protein-ligand docking, motifscore, was developed. it is non-energy-based, and docking is, instead, scored by counting the occurrences of motifs of protein-ligand interaction networks constructed using structures of protein-ligand complexes. motifscore has been tested on a benchmark set established by others to assess its ability to identify near-native complex conformations among a set of decoys. in this benchmark test, 84% of the highest-scored docking conformations had root-mean-square deviations below Å from the native conformation, which is comparable with the best of several energy-based docking scoring functions. many of the top motifs, which comprise a multitude of chemical groups that interact simultaneously and make a highly significant contribution to motifscore, capture recurrent interacting patterns beyond pairwise interactions.  while providing quite good docking scores, motifscore is quite different from conventional energy-based functions. motifscore thus represents a new, network-based approach for exploring problems associated with molecular docking. 