the canonical glutathione transferases ( also known as gsts ; ec 
 2.5.1.18 ) catalyze addition of an excellent nucleophile to an electrophilic 
 center .
they play important roles in the metabolism and detoxification 
 of numerous endogenous and xenobiotic compounds , including oxidized 
 lipids , drugs , and pollutants .
they consist of an n - terminal 
 thioredoxin domain and a c - terminal -helical domain .
although 
 a number of gsts are known to have specific substrates , many if not 
 most have no clearly assigned biological substrates .
in addition , 
 the general nature of their chemistry leads to enzymes that tend to 
 be catalytically promiscuous even with respect to the transition state 
 for the reaction . as a consequence ,
computational docking methods have been 
 widely used in ligand discovery 
 for many enzymes . in particular
, these methods can predict the docking pose of a known 
 ligand ( docking ) and/or predict ligands in a large library of small 
 molecules ( virtual screening ) . docking consists of searching through 
 plausible binding modes of a compound and scoring each mode to distinguish 
 a near - native binding pose from others ( docking ) .
virtual screening 
 consists of performing docking for each candidate ligand , followed 
 by the ranking of the candidate ligand by their best docking scores .
although this structure - based approach has been used successfully 
 for the prediction of both substrates and other types of ligands ( e.g. , 
 orthosteric inhibitors and allosteric modulators ) , docking for substrate 
 discovery is most difficult , particularly when the enzyme may accommodate 
 more than one type of reaction .
in addition , only ground state or 
 intermediate state complexes , not transition states , are typically 
 accessible by standard docking procedures .
the resulting complexes 
 may or may not be directly competent for turnover in the absence of 
 information on the  preorganization  of the enzyme 
substrate 
 complex required for catalysis . despite the difficulties involved 
 in the docking of substrates ,
structure - based methods have been used 
 successfully for substrate discoveries in several systems .
there are two principal challenges in the de novo prediction of the substrate preferences of enzymes in large , functionally 
 diverse superfamilies .
the first challenge is the availability of 
 experimentally determined enzyme structures , which lags behind the 
 number of known protein sequences by a factor of approximately 400 
 as of june 2013 .
to remedy this situation , virtual screening has also 
 relied on comparative models , not only experimentally determined structures .
the second challenge is that standard docking procedures lack sufficient 
 constraints that efficiently define the productive geometries between 
 the reacting species ( substrates ) on the enzyme surface .
this issue 
 is addressed here by applying a covalent bond constraint between the 
 sulfur of gsh and the electrophilic substrate .
the 
 absence of a large - scale benchmarking of covalent docking using 
 either x - ray structures or comparative models , and the need to predict 
 substrates for gst enzymes , inspired us to investigate the following 
 questions .
can covalent docking accurately predict docking poses of 
 known ligands , given experimentally determined , homologous , or modeled 
 structures in the gst superfamily ? can covalent docking accurately 
 predict ligands despite the catalytic promiscuity of many gst enzymes ? 

what is the difference in the utility of apo , holo , comparative modeling , 
 and homologous structures for virtual screening in the gst superfamily ? 
 can the virtual screening be improved by consensus scoring , relying 
 on independent screening against multiple holo , apo , comparative modeling , 
 and homologous structures in the gst superfamily ?
if multiple models 
 are calculated on the basis of different templates , can any of them 
 outperform apo and even holo x - ray structures of the target ?
if so , 
 can one reliably identify which model will do so , or even perform 
 optimally among a set of modeled structures ; are there sequence and/or 
 structural attributes ( i.e. , the overall target  template sequence 
 identity , the binding site target  template sequence identity , 
 and the predicted accuracy of a model ) that reliably predict the accuracy 
 of ligand docking ? in this report , we attempt to answer these 
 questions with the aid 
 of a virtual library of compounds and 14 representative gst enzymes 
 of known structure and function .
the virtual library consists of known 
 substrates for the selected gst proteins , and a large number of compounds 
 selected from kegg and metacyc .
for each target , the entire virtual library 
 of compounds was docked to the known x - ray structures , homologous 
 structures , and comparative models of the protein .
the results are 
 analyzed by comparing the docking poses of native products to x - ray 
 structures and calculating the enrichment of known products with respect 
 to the entire compound library .
we begin by describing the gst 
 catalyzed reactions , the docking 
 library , the selected gst proteins , the automated modeling pipeline , 
 the docking pipeline , and methods to evaluate the accuracy of predicted 
 docking poses for known ligands and the accuracy of virtual ligand 
 screening ( methods ) .
we then describe and 
 compare the results of docking native ligands and virtual screening 
 using apo , holo , comparative modeling , and homologous structures ( results ) .
finally , we discuss the implications of 
 the current approach and answer the questions we asked previously , 
 given our modeling , docking , and benchmark ( discussion ) .
we begin by listing the set of the gst catalyzed 
 reactions used 
 in this study , followed by a description of the corresponding products 
 that comprise a virtual screening library .
next , we describe how we 
 selected gst targets for docking and how we built their comparative 
 models .
finally , we describe the covalent docking pipeline as well 
 as the assessment criteria for evaluating the accuracy of docking 
 and virtual screening .
, we considered 11 different 
 kinds of gst catalyzed reactions , each with a different substrate 
 motif ( figure 1 ) . for each substrate motif , 
 we found all matching molecules in the kegg and metacyc compound databases . 

the search was performed with the oechem tools by matching a smiles string 
 of a database molecule to the smarts strings 
 representing the substrate motifs .
we also added 64 substrates from the literature , for the total of 6,738 unique substrates .
each row 
 defines a type of reaction based on the reactive functional group . 

next , each of the substrates was converted to one or more 
 products , 
 as follows .
the conversions corresponding to the 11 reactions were 
 carried out using oechem s library generation function with explicit hydrogens .
the reactive functional 
 groups in the substrate were identified by comparing the substrate 
 smiles string with the smarts string representing the substrate motif 
 undergoing conversion to a product motif , for each reaction .
each 
 match was then used to convert the substrate to a product that was 
 added to the virtual screening library . finally , we prepared 
 products for docking . for products with undefined 
 stereocenters ,
we first enumerated the stereoisomers using oechem , with the maximum number of stereoisomers retained 
 arbitrarily set to 16 for computational efficiency .
second , protonation 
 states of each stereoisomer was enumerated within ph range 68 , 
 followed by generating a conformation for each protonation state , 
 using epik and ligprep .
finally , force field parameters for each product were generated 
 using the hetgrp_ffgen utility ( schrdinger , lcc ) .
a subset of gst structures 
 was selected from the protein data bank ( pdb ) for retrospective docking 
 by maximizing the sequence and functional coverage of the gst superfamily 
 ( figure 2 ) , following three steps : first , all 
 gst structures with a cocrystallized gsh - substrate conjugate were 
 extracted from the pdb .
second , structures for the same protein were 
 grouped together , resulting in 14 groups with unique sets of ligands . 
 finally , for each group , the structure with the highest resolution , 
 the lowest rfree , , and a unique ligand 
 was selected , resulting in 14 target holo x - ray structures ( supporting information table s1 ) .
representative network 
 view of the gst superfamily showing where 
 docking targets fall in this superfamily .
in this view , 2,190 50% 
 sequence identity filtered nodes (  id50 node  ) representing 
 13,493 sequences are shown as dots , and each sequence similarity relationship 
 with a blast e - value  1  10 is shown as a line or edge between representative nodes .
the 13,493 sequences 
 for cytosolic gsts (  cytgsts  ) are based on pfam ( v26.0 ) sequences that were at least 100 residues in 
 length and had scores above the gathering threshold for at least one 
 pfam hidden markov model ( hmm ) corresponding to the cytgst n - terminal 
 domain (  gst_n ,  
also included in the data set were 58 
 proteins that lacked matches to a gst_n * hmm but were chosen as enzyme 
 function initiative ( efi ) targets for 
 experimental characterization because of their sequence similarities 
 to cytgsts . sequence similarity was detected with an all - by - all blast 
 search using blastp ( v2.2.24 ) .
sequence identities are calculated 
 using cd - hit (  id50  set ) . 

the edges are shown using the organic layout where nodes that are 
 more highly interconnected are clustered more tightly together .
pdb 
 structures associated with cytgst sequences ( 95% or more sequence 
 identity to the pdb sequence ) were identified using the sfld interface .
an id50 node is shown as a blue dot 
 if any member sequence of the node was associated with a pdb structure . 

representative nodes where only an efi structure was associated with 
 its member sequences are shown as red dots .
retrospective docking 
 targets are marked using dark blue borders ; only 10 targets are shown 
 because the rest are more than 50% identical to these 10 representatives 
 shown . because gsts function as dimers , 
 we used for docking the biological 
 dimer unit from the pdb .
the substrate part of the gsh - substrate conjugate 
 in holo structures was removed , resulting in gst holo x - ray structures 
 with only gsh bound .
each dimer was then prepared for docking using 
 the protein preparation wizard ( schrdinger , lcc ) by designating 
 each of the two gshs as the reactive cofactor , resulting in two receptor 
 binding sites ( with the exception of 1gwc , which has only one site 
 occupied by gsh ) ; in most cases , the two binding sites are slightly 
 different because symmetry is generally not imposed during crystallographic 
 structure refinement .
comparative 
 models for each target were built by selecting template dimer structures 
 from the pdb with either an unmodified gsh ( apo ) or a conjugate form 
 ( product ) of gsh ( holo ) .
for each target , we aimed to include one 
 holo and one apo template structure for each of the four sequence 
 identity ranges : 2540% , 4055% , 5570% , and 
 7085% .
the target  template sequence alignments were 
 generated using hhalign , which in turn 
 relied on target and template profile hidden markov models ( hmms ) 
 generated using hhblits with the uniprot20 
 database .
five - hundred models were built 
 for each target  template alignment , using the automodel protocol 
 in modeller-9v10 .
the coordinates for 
 gsh atoms in the template structure were then copied to the target 
 models , using the blk function of modeller .
models were then assessed 
 by the discrete optimized protein energy ( dope ) function of modeller .
finally , the model with the lowest 
 normalized dope score was used for docking . to compare different 
 comparative models and compare the models against x - ray structures , 
 we calculated the binding site sequence identity of the templates 
 to the x - ray structures , and the binding site non - hydrogen atom root - mean - square 
 deviation ( rmsd ) between a comparative model and x - ray structure . 

the binding site was defined to contain atoms in all residues with 
 at least one atom within 6  of any ligand atom in the x - ray 
 structure .
as expected , strong correlations were observed between 
 the overall target  template sequence identity and the target 
template 
 binding site sequence identity and the model s non - hydrogen 
 atom rmsd error ( figure 4a , b ) .
as designed by the choice of the template 
 structures , the distribution 
 of the non - hydrogen atom rmsd error of the binding site confirms that 
 there is a range of accuracy across the set of modeled structures . 

the models with subangstrom non - hydrogen atom rmsd error tend to be 
 those based on templates with over 80% sequence identity in the binding 
 site residues to the target sequence . over a third of models had a 
 binding site rmsd error over 4  .
thus , homology models cover 
 a range of sequence identities and a range of accuracy .
covalent docking of a 
 potential product to a receptor was started by placing it in the binding 
 site .
more specifically , the coordinates of the gsh part of the potential 
 product were matched to the gsh molecule in the receptor , and the 
 coordinates of the remaining atoms of the product were then built 
 using omega .
up to 20 initial configurations 
 were generated for each product molecule . for each initial configuration , 
 we then used plop s tether pred ( academic version 25.6 ) function to rotate the rotatable bonds in the gsh - substrate 
 conjugate
we sampled 
 all of the rotatable bonds in the substrate part of the gsh - substrate 
 conjugate as well as the ca  cb and cb 
sg bonds in the 
 cysteine residue of gsh . up to 50 configurations were generated by 
 plop starting from each initial configuration .
these configurations 
 were then scored by calculating the total potential energy of the 
 product  receptor complex by plop ( in the units of kilocalories 
 per mole ) , which in turn relies on the opls force field with a variable - dielectric 
 generalized born model . the product 
receptor 
 distances alone were also scored by the atomic statistical potential 
 posescore ( in arbitrary units ) . while 
 the plop potential energy does not contain the entropic contributions 
 to the binding free energy , the posescore term does approximate the 
 contribution of the interface to the binding free energy .
finally , for virtual screening , a substrate 
 was ranked using the 
 median plop energy of its products different stereoisomers 
 and protonation states .
the plop energy of a product in a specific 
 stereoisomer form and protonation state was the median plop energy 
 of its different configurations . when a 
 native product was docked ( supporting information table s2 ) , the docking pose of the product
was assessed for accuracy 
 based on its non - hydrogen atom rmsd from the native configuration , 
 after superposition of the receptor used for docking on the native 
 structure of the receptor ; gsh was excluded , except for the cysteine 
 sulfur atom .
the accuracy of virtual screening was evaluated 
 by the enrichment for the known products ( supporting 
 information table s3 ) among the top scoring potential products . 

the enrichment curve was obtained by plotting the percentage of actual 
 products found ( y - axis ) within the top ranked subset 
 of all database compounds ( x - axis on logarithmic 
 scale ) .
logauc , the area under the curve of the enrichment plot , was 
 also calculated to indicate the accuracy of enrichment ; random selection 
 has a logauc of 14.5 .
we begin by evaluating the 
 docking pose of the native ligands using 
 holo and apo structures , followed by evaluating the docking pose of 
 the native ligands using comparative models and homologous structures . 

next , enrichment of the known ligands using x - ray structures is benchmarked 
 and analyzed .
finally , we analyze the enrichment of the known ligands 
 using comparative models . as the easiest 
 test
, we first applied covalent docking by plop to reconstruct binding 
 poses of the crystallographic products in the corresponding holo x - ray 
 structures ( supporting information table 
 s1 ) . in all 14 cases , a docking pose within 3  of the native 
 structure was sampled ( within 2  for 13 cases ) , though not necessarily 
 recognized as such by the scoring function ( supporting 
 information table s4 ) ; the 3  threshold on the all - atom 
 rmsd subjectively distinguishes near - native and nonnative poses . however , 
 the top 3 ranked structure was within 3  to the native structure 
 in 9 out of 14 cases using the plop energy function , and in 11 out 
 of 14 cases using posescore ( figure 3 ) .
next , we performed a slightly more difficult , but still easy , test . 
 for targets with more than one crystallographic product
, we docked 
 each product to the nonnative holo x - ray structure of the target ( cross - docking ; supporting information table s5 ) . in all 12 cases , 
 a docking pose within 3  to the native structure was sampled 
 ( within 2  for 11 cases ) .
the top 3 ranked complex was within 
 3  to the native structure in 7 out of 12 cases for the plop 
 energy function , and in 8 out of 12 cases for posescore . to test
whether 
 combining the two scoring functions can further improve the result , 
 we considered an optimal linear combination of the posescore and the 
 plop potential energy , using the average all - atom rmsd of the top 
 ranked docking poses as the optimization criterion ; the optimal linear 
 combination has zero weight for the plop potential energy . finally , 
 to test whether or not apo structures can be used for 
 predicting the native product docking pose , we docked the two products 
 of gstm1_human ( gdn and gtd ) to its two available apo structures ( 1xw6 and 1yj6 ) .
the results were 
 comparable to those from cross - docking ( supporting 
 information table s6 ) .
comparison of the predicted complex between 
 the native product 
 and its holo x - ray structure with the native pose : ( a ) 1f3b with gbx ; ( b ) 2f3 m with gtd ; ( c ) 2gst with gps ; ( d ) 1m9b with ibg .
the surface 
 of the receptor is shown , with red and blue corresponding to oxygen 
 and nitrogen atoms .
the product is shown in the native configuration 
 ( gray ) , the most accurate sampled configuration ( yellow ) , and the 
 top ranked configuration by the plop energy function ( light blue ) 
 and posescore ( green ) . to test whether or not comparative 
 models can be used for docking pose prediction of native products , 
 a total of 62 homology models were generated for the 14 gst targets 
 ( supporting information table s7 ) .
we performed 
 native product docking against these models as well as some of their 
 template structures . to find out to what degree
the success 
 of docking pose prediction was limited by sampling versus scoring , 
 we plotted the rmsd error of the top ranked docking pose versus the 
 rmsd error of the best sampled docking pose ( figure 4c ) .
as expected there 
 was a significant correlation between these two rmsds ; the rmsd error 
 of the best sampled pose is a lower bound of the rmsd error of the 
 top ranked pose .
even when accurate poses were sampled , the scoring 
 function often did not recognize them ( figure 4c ) . in 56 out of 62 cases , a docking pose within 3  was sampled , 
 but only in 22 out of 62 cases , it was ranked in the top 3 .
therefore , 
 the limiting factor for docking pose prediction is the accuracy of 
 the scoring function used . in this application
, the inaccuracy of 
 the scoring function may be exacerbated by errors in the model of 
 the binding site ( figure 6 ) . despite the difficulties 
 in ranking a near - native pose in the context of comparative model , 
 for 6 of the 14 targets ( 1gwc , 2ab6 , 3ljr , 3o76 , 2c4j , and 1bye ; figure 5 ) , 
 the rmsd errors of the top ranked poses by posescore were comparable 
 to those against holo x - ray structures .
we then asked the question 
 whether or not there were model properties 
 that predict the best model for docking pose prediction , among multiple 
 models based on different template structures .
there was essentially no correlation between a model s binding 
 site rmsd error and the rmsd error of the top ranked docking pose 
 ( figure 4d ) .
there were only weak correlations 
 between the target  template sequence identity and the rmsd 
 error of the best - sampled docking pose ( figure 4f ; r = 0.19 for apo tempalates , and r = 0.20 for holo templates ) .
we also calculated 
 the correlation between the rmsd errors of the model binding site 
 and the best - sampled docking pose .
we also investigated whether or not 
 holo and apo structures performed differently for docking pose prediction . 
 on average , models based on holo template structures were not noticeably 
 better for docking pose prediction than models from apo template structures .
finally , we tested the accuracy of docking pose prediction using 
 a homologous structure itself , instead of a comparative target model 
 based on it .
the test was performed with the two known products of 
 gstm1_human ( gdn and gtd ) and its homologous structure ( 2c4j ) .
the docking pose 
 of gtd had the rmsd error of 2.2  , and the docking pose of gdn 
 had the rmsd error of 4.7  , on average not significantly different 
 from that using the corresponding comparative model ( 3.0  for 
 gtd and 3.5  for gdn , binding site rmsd error ) .
comparison of comparative 
 model properties and the rmsd errors 
 of the docked poses .
each point represents the property of a single 
 comparative model and/or the rmsd error of the docked pose of its 
 corresponding crystallographic product against this model .
results 
 for models built based on apo and holo template structures are shown 
 in red and blue , respectively , with the least - squares linear fits 
 shown as dashed lines .
the gray dashed line indicates the lower 
 bound of the rmsd of the top 3 ranked poses .
non - hydrogen atom rmsd error of the docked poses , for the crystallographic 
 product against its corresponding comparative models with different 
 sequence identities . for each target , the rmsd errors of the best - sampled 
 poses using apo templates ( stripped red bar ) , the rmsd errors of the 
 top - ranked poses by posescore using apo templates ( stripped blue bar )
, 
 the rmsd errors of the best - sampled poses using holo templates ( red 
 bar ) , and the rmsd errors of the top ranked pose by posescore using 
 apo templates ( blue bar ) are shown .
for comparison , the rmsd errors 
 of the docked poses using holo x - ray structures are shown as horizontal 
 lines ( blue solid line for the top ranked pose by posescore , red dashed 
 line for the best - sampled pose ) .
( a ) successful 
 docking of gbx onto a model of 3csh , based on 85% sequence identity 
 to the template 3o76 ( chain a ) .
the ribbon representation of the crystal structure is 
 shown in light gray and that for the model in orange .
the c- 
 rmsd error of the model is 0.64  , and the non - hydrogen atom 
 rmsd error of the binding site is 0.71  .
the native pose is shown in light gray , 
 the best - sampled pose ( rmsd of 0.90  ) is shown in yellow , the 
 pose top ranked by plop ( 0.92  ) is shown in cyan , and the pose 
 top ranked by posescore ( 0.92  ) is shown in green ( not visible 
 because it is hidden behind the cyan structure ) .
( b ) failed docking 
 of gps onto a model of 1f3b , based on 76% sequence identity to the template 1ydk ( chain b ) .
the ribbon 
 representation of the crystal structure is shown in light gray and 
 that for the model in orange .
the c- rmsd error of the model 
 is 2.42  , and the non - hydrogen atom rmsd error of the binding 
 site is 2.47  .
the native pose is shown in light gray , the best - sampled pose ( rmsd 
 of 5.38  ) is shown in yellow , the pose top ranked by plop ( 6.58 
  ) is shown in cyan , and the pose top ranked by posescore ( 0.92 
  ) is shown in green .
we now address the question whether or not holo 
 and apo x - ray structures can be used for virtual product screening 
 by covalent docking .
we first applied plop to virtual substrate 
 discovery by docking the whole product library against holo structures 
 of 11 gsts with more than 3 known substrates ; this is the minimal 
 number of known ligands needed for computing an enrichment curve . 

we docked the products onto each chain of the gst dimers independently 
 and analyzed the docking results separately for each chain as well 
 as collectively for all chains , by relying on our consensus scoring 
 scheme ; consensus scoring combines docking 
 scores from independent virtual screen against multiple chains / structures / models 
 of the same target . in all 11 cases ,
the enrichment result was better 
 than that from random selection ( figure 7 ) . 
 in 7 of the 11 cases , logauc was better than 25 .
in contrast to the 
 noncovalent docking benchmark , consensus 
 scoring did not improve logauc for any target .
next , we tested 
 the enrichment of known ligands using two apo structures 
 of gstm1_human ( 1yj6 and 1xw6 ; supporting information figure s2 ) .
the averaging 
 logauc ( 32.6 ) is slightly worse than that using the corresponding 
 holo structure ( 33.3 ) . finally , we tested whether or not consensus 
 scoring using multiple 
 structures could improve enrichment for known products by combining 
 virtual screening results from either two apo or two holo structures 
 of gstm1_human . using the two apo structures ( 1yj6 and 1xw6 ) ,
similarly , using holo structures ( 2f3 m and 1xwk ) resulted in a logauc 
 value of 33.6 . in both cases , logauc for consensus scoring was between 
 the best and worst logauc values of individual virtual screening runs .
enrichment 
 curves for the 11 gst targets with 3 or more known products , 
 using holo x - ray structures . for each target , the enrichment curves 
 for the individual chains ( blue , green , yellow , and cyan lines ) , the 
 consensus scoring for multiple chains ( red line ) , and random selection 
 ( black dashed line ) are shown . to evaluate 
 the utility of homology models for discovering substrates of gsts , 

the enrichment of known products was determined and compared to that 
 for the crystal structures .
we first tested whether or not there 
 was an apparent correlation between the enrichment and the target 
template 
 sequence identity by virtual screening against the 6 comparative models 
 of the target 2f3 m .
the absence of 
 such a correlation was also observed in an extensive noncovalent docking 
 benchmark .
given the absence of 
 a correlation between the enrichment and the 
 target  template sequence identity , we next benchmarked the 
 utility of the best comparative models for virtual screening , for 
 the 11 targets with more than 3 known substrates ; the best comparative 
 model is that based on the highest target  template sequence 
 identity .
we docked the products onto each chain of the gst dimers 
 independently , and analyzed the docking results separately for each 
 chain as well as collectively for all chains , by relying on our consensus 
 scoring scheme .
in all cases , the enrichment 
 result was better than that from random selection ( figure 9 ) . in 5 of the 11 cases , logauc was better than 
 25 .
we further tested whether or not consensus docking using multiple 
 comparative models can improve the enrichment by combining the docking 
 results from the 6 comparative models for target 2f3 m ( figure 8) .
consensus scoring resulted in a logauc of 40.4 , 
 which was slightly better than the best logauc value for individual 
 models ( 40.3 ) and was better than logauc for the holo x - ray structure .
finally , we tested whether or not homologous x - ray structures could 
 be used in place of comparative models for virtual screening , using 
 a homologue of gstm1_human ( 2c4j ) .
the corresponding logauc is 30.9 
 for chain a , 31.6 for chain b , and 31.5 for consensus scoring using 
 chains a and b , which is significantly worse than that for the comparative 
 model based on 2c4j ( 34.9 for chain a , 35.8 for chain b , and 38.7 for consensus scoring ) . 

enrichment curves for target 2f3 m using its 6 comparative models based 
 on different apo and holo templates . for each model ,
the enrichment 
 curves for the individual chains ( blue , green , yellow , and cyan lines ) , 
 the consensus scoring for multiple chains ( red line ) , and random selection 
 ( black dashed line ) are shown .
template sequence identities 
 are as follows : the 2f3m-3t2u and 2f3m-2fhe sequence identities 
 are 33% and 50% , respectively ( apo templates ) .
template 
 sequence identities are as follows : 33% ( 2f3m-3t2u , apo ) , 50% ( 2f3m-2fhe ) , 33% ( 2f3m-2gsr , holo ) , 
 42% ( 2f3m-1m99 , holo ) , 66% ( 2f3m-1gsu , holo ) , and 85% 
 ( 2f3m-2c4j , holo ) .
binding 
 site rmsd errors are as follows : 5.9  ( 2f3m-3t2u , apo ) , 3.1  
 ( 2f3m-2fhe ) , 7.4  ( 2f3m-2gsr , holo ) , 3.1  
 ( 2f3m-1m99 , holo ) , 1.4  
 ( 2f3m-1gsu , holo ) , and 1.1 
  ( 2f3m-2c4j , holo ) .
enrichment curves for the 11 gst targets with 
 3 or more known products , 
 using their comparative models based on templates with the highest 
 sequence identities . for each target , the enrichment curves for the 
 individual chains ( blue , green , yellow , and cyan lines ) , the consensus 
 scoring for multiple chains ( red line ) , and the random selection ( black 
 dashed line ) are shown .
template sequence identities 
 are as follows : 85% ( 2f3m-2c4j ) , 85% 
 ( 3csh-2c4j ) , 77% ( 2c4j-6gsv ) , 54% ( 3ljr-2c3q ) , 33% ( 2gsq-2gsr ) , 28% ( 3gx0 - 1pn9 ) , 76% ( 2gst-2c4j ) , 44% ( 1m9b-2c4j ) , 64% ( 1vf3 - 3ktl ) , 76% ( 2ab6 - 6gsv ) , and 46% ( 1gwc-2vo4 ) .
we developed a comparative 
 modeling and covalent docking pipeline 
 for predicting docking poses of known products and for the virtual 
 discovery of substrates , as well as benchmarked it on the gst enzymes . 

either a holo or apo structure can be used for predicting docking 
 poses relatively accurately in a majority of cases .
while using homology 
 models for predicting docking poses is more difficult , a homology 
 model can also be useful . specifically , for native holo x - ray structures 
 as receptors , one of the top 3 ranked models is within 3  all - atom 
 rmsd of the native complex in 11 of the 14 test cases ( in 8 of 12 
 cases for cross - docking in which a product from another holo structure 
 of the same enzyme is docked ; supporting information tables s4 and s5 ) .
for comparative models based on more than 30% 
 sequence identity , one of the top 3 ranked models is within 3  
 all - atom rmsd of the native complex in 22 of the 62 test cases ( figure 5 ) . for virtual screening against holo x - ray structures , 

the enrichment logauc value is better than random ( 14.5 ) in all cases 
 and better than 25 in 7 of 11 cases ( figure 7 ) . for models based on templates with the highest sequence identity 
 ( approximately 3085% ) ,
the logauc is better than 25 in 5 of 
 11 cases ( figure 9 ) , not significantly different 
 from the crystal structures . in conclusion , we show that covalent 
 docking can be a useful tool for substrate discovery .
first , we assess to what 
 degree the accuracy of the docking poses and ranking of alternative 
 products is limited by the degree of structure sampling and the accuracy 
 of scoring .
second , we discuss the utility of comparative models for 
 docking and virtual screening .
third , we conclude by commenting on 
 docking transition state conformations of ligands . to predict an accurate docking 
 pose ,
two conditions need to be satisfied : ( i ) a near - native structure 
 needs to be generated by the sampling procedure and ( ii ) an accurate 
 sample structure needs to be scored better than the inaccurate structures 
 by the scoring function . using x - ray structures as receptors ,
a near - native 
 docking pose ( <3  ) was sampled in all cases ( supporting information tables s4 and s5 ) ; thus , 
 the degree of sampling does not limit the prediction accuracy .
but , 
 in 21% of the native docking cases and in 33% of the cross - docking 
 cases , a near - native model was not ranked in top 3 by posescore . using 
 homology models as receptors , a near - native docking pose
was sampled 
 in 56 of 62 cases but was only ranked in top 3 in 22 of 62 cases ( figure 5 ) .
the reason is that our scoring function is not 
 always accurate enough to differentiate between accurate and inaccurate 
 docking poses . to rank different substrates
we chose the 
 median instead of the lowest plop energy because it resulted in a 
 significantly higher average logauc ( 27.3 vs 18.9 ) .
a likely reason 
 is that the energy estimate of a single configuration is noisier than 
 the median , compensating for the lack of accounting for the proper 
 physics of the problem . in the future
, other functions such as the 
 boltzmann average of the sampled energies will be explored .
although posescore reranking of the docking poses generated by 
 optimizing the plop energy further improved the accuracy of the top 
 scoring docked poses , this result does not imply that posescore is 
 more accurate than the plop energy .
it is possible that reranking 
 by plop would also improve the accuracy of the poses generated by 
 optimizing posescore ; for technical reasons , it is not straightforward 
 to actually make this test .
our result indicates merely that the two 
 scoring functions are different and imperfect and that posescore contains 
 helpful information not present in plop , while being silent on whether 
 or not plop also contains information not present in posescore .
in contrast to noncovalent docking , on 
 average , consensus scoring did not improve virtual screening 
 accuracy .
a combination of the following two assumptions explains 
 this finding : ( i ) covalent docking is limited more by scoring than 
 by sampling compared to noncovalent docking ; and ( ii ) consensus scoring 
 works best when sampling is a limiting factor in docking .
a major 
 difference between noncovalent and covalent docking is that the latter 
 is constrained by a covalent bond , resulting in an easier sampling 
 problem all other things considered equal .
using multiple receptor 
 structures to mimic receptor flexibility was already shown to lead 
 to more accurate sampled docking poses in noncovalent docking .
thus , the disadvantage of a larger number of 
 decoys produced by using multiple templates , which increases the burden 
 on the scoring function for identifying the most accurate pose , is 
 outweighed by the more accurately sampled poses that are more likely 
 to be ranked highly by the scoring function .
in contrast , if scoring 
 is limiting and sampling is not , as is the case in covalent docking 
 ( cf . our native docking and cross - docking tests ) , considering poses 
 from multiple receptor structures is not likely to improve the results , 
 as observed . 
given the demonstrated limitations of current 
 scoring , the accuracy 
 of ranking ligand poses and ligands can be improved by more accurate 
 scoring functions . on the one hand , the ability of the plop potential 
 energy to rank the native state ( or at least a near - native state ) 
 may be improved by estimating the free energy of the native state .
on the other hand , the accuracy of the statistical 
 potential could be improved by applying the bayesian framework for 
 inferring statistically optimized atomic potentials ( soap ) .
this framework ( 1 ) uses multiple data - driven 
  recovery  functions based on probability theory , without 
 recourse to questionable statistical mechanical assumptions about 
 statistical potentials ; ( 2 ) restrains the relative orientation between 
 two covalent bonds instead of a simple distance between two atoms , 
 in an effort to capture orientation - dependent interactions such as 
 hydrogen bonds ; ( 3 ) performs bayesian smoothing for estimating the 
 underlying smooth distributions from noisy observations ; ( 4 ) applies 
 bayesian inference for calculating parameter values that maximize 
 the posteriori probability ; and ( 5 ) benefits from bayesian model selection 
 based on bayesian predictive densities , in an effort to improve the 
 generalizability of the derived statistical potentials . using comparative 
 models for predicting docking poses
one possible reason is that covalent docking is sensitive 
 to minor changes in the binding site , especially for residues close 
 to the gsh sulfur atom due to the substrate 
another possible 
 reason is that the active sites of gst enzymes are often defined in 
 part by the c - terminus , which is often flexible or missing in a template 
 structure ( figure 6 ) .
moreover , there was no 
 significant correlation between the binding site non - hydrogen atom 
 rmsd and the docking pose rmsd ( figure 4 ) . 

the lack of this correlation may also result from the disproportionate 
 impact of residues close to the gsh sulfur atom on the docking pose .
virtual screening using comparative models generated results comparable 
 to the x - ray structures ( average logauc of 28.2 versus 27.3 , better 
 than 25 in 5 versus 7 out of 11 cases ) , consistent with what has been 
 observed when using comparative models for virtual screening by noncovalent 
 docking ( average logauc of 28.7 versus 30.6 ) , even though comparative models performed much worse than x - ray structures 
 for docking pose prediction ( average non - hydrogen atom rmsd of 3.98 
  versus 1.83  ) .
this finding may be rationalized by comparatively 
 worse scoring of both true and decoy products ( average posescore of 
 4.2 and 4.1 versus 60.0 and 52.2 ) when docking 
 against inaccurate binding sites in comparative models . in prospective 
 applications , if some substrates for the target 
 of interest are known , virtual screening using comparative models 
 can be improved by selecting the most enriching model for known substrates . if no substrates are known for the target of interest , it might 
 be advantageous to build and dock into multiple structurally diverse 
 models and analyze the virtual screening results by hand .
an experienced user may be able to identify 
 an accurate docking result based on the totality of their knowledge 
 about the proteins of interest .
generally , a substrate of an enzyme 
 is required to have a sufficiently 
 high kcat / km ; a typical threshold is 10 ( mol / l ) s. an enzymatic reaction involves interconversions 
 between pairs of stable states , separated by a transition state ; the 
 stable states include the enzyme and substrate without interaction 
 ( e + s ) , an enzyme  substrate complex intermediate ( es ) , an 
 enzyme  product intermediate ( ep ) , and the enzyme and product 
 without interaction ( e + p ) . in general , kcat
thus , virtual ligand screening should 
 in principle consider all states , which is generally not the case .
when conversion from es to ep is the rate - limiting step , kcat is mostly determined by the transition state 
 energy barrier for this step .
enzymes lower this barrier by creating 
 a binding site that is complementary to the transition state of the 
 substrate .
thus , docking a substrate ( or a product ) in its transition 
 state conformation may be more appropriate than docking a substrate 
 ( or a product ) in its ground state .
although quantum mechanics calculations 
 can model the transition state and estimate the activation energy 
 barrier , running these calculations for the entire virtual screening 
 library is not computationally feasible .
however , the enrichment results could possibly be 
 greatly improved if we can dock the transition state , or a transition 
 state like structure , without recourse to quantum mechanics . for gst 
 substrate discovery ,
compounds similar to the high - energy intermediates 
 used for docking to amidohydrolases could 
 be constructed , provided some technical challenges are solved ( e.g. , 
 building c atoms with 2 partial bonds for a total of 5 bonds ) .
although 
 enzymes presumably maximize complementarity to the transition state , 
 other states during the reaction must also be compatible with the 
 potentially flexible binding site of the enzyme .
thus , substrate discovery 
 might also be improved by consideration of docking of all ligand states 
 to the corresponding state(s ) of the enzyme .
can covalent docking accurately predict docking poses of 
 known ligands , given experimentally determined , homologous or modeled 
 structures in the gst superfamily ?
yes , either a holo or apo 
 structure can be used for predicting docking poses accurately in a 
 majority of cases , while using homology models is less successful . 

specifically , for native holo x - ray structures as receptors , one of 
 the top 3 ranked models is within 3  all - atom rmsd of the native 
 complex in 11 of the 14 cases ( in 8 of 12 cases for cross - docking ) . 
 for comparative models as receptors ,
one of the top 3 ranked models 
 is within 3  all - atom rmsd of the native complex in 22 of the 
 62 test cases .
can covalent docking accurately predict 
 ligands despite the 
 catalytically promiscuity of many gst enzymes ? often .
what is the difference in the utility of apo , holo , comparative 
 modeling , and homologous structures for virtual screening in the gst 
 superfamily ?
holo and apo structures , homologous structures , 
 and comparative models can all enrich for known products . for virtual 
 screening against holo x - ray structures ,
the enrichment logauc value 
 is better than random in all cases , and better than 25 in 7 of 11 
 cases ( average logauc 27.3 ; figure 7 ) .
using 
 apo x - ray structures of 2f3 m , the average logauc ( 32.6 ) is slightly worse than 
 that using the corresponding holo structure ( 33.3 ) . for models based 
 on templates with the highest sequence identity ( approximately 3085% ) , 

the logauc is better than 25 in 5 of 11 cases ( average logauc 28.2 ; 
 figure 9 ) , which is comparable to that for 
 x - ray structure .
using the homologous structure of 2f3 m , the logauc ( 31.3 ) 
 is significantly worse than that using comparative models ( 36.5 ) .
can the virtual screening be improved by consensus scoring , 
 relying on independent screening against multiple holo , apo , comparative 
 modeling , and homologous structures in the gst superfamily ? rarely .
consensus scoring using multiple chains from the same structure / model 
 or multiple structures / models of the same target often does not improve 
 the virtual screening accuracy .
although consensus scoring improves 
 logauc in comparison to using a single receptor in some cases , it 
 makes it worse in others .
if multiple models are calculated 
 on the basis of different 
 templates , can any of them outperform apo and even holo x - ray structures 
 of the target ?
if so , can one reliably identify which model will do 
 so , or even perform optimally among a set of modeled structures ; are 
 there sequence and/or structural attributes ( i.e. , the overall target  template 
 sequence identity , the binding site target  template sequence 
 identity , and the predicted accuracy of a model ) that reliably predict 
 the accuracy of ligand docking ? for docking pose prediction , 
 there are a few cases where comparative models outperformed the corresponding 
 holo x - ray structure ( figure 5 ) .
however , we 
 did not find a predictor that could reliably predict the accuracy 
 of the docked pose . for virtual screening , some models also outperform 
 the corresponding holo x - ray structure ( figures 8 and 9 ) , but we did not find a correlation 
 between logauc and the binding site sequence identity or the binding 
 site rmsd errors .
 enzymes 
 in the glutathione transferase ( gst ) superfamily catalyze 
 the conjugation of glutathione ( gsh ) to electrophilic substrates . 
 as a consequence 
 they are involved in a number of key biological processes , 
 including protection of cells against chemical damage , steroid and 
 prostaglandin biosynthesis , tyrosine catabolism , and cell apoptosis . 
 
 although virtual screening has been used widely to discover substrates 
 by docking potential noncovalent ligands into active site clefts of 
 enzymes , docking has been rarely constrained by a covalent bond between 
 the enzyme and ligand . in this study 
 , we investigate the accuracy 
 of docking poses and substrate discovery in the gst superfamily , by 
 docking 6738 potential ligands from the kegg and metacyc compound 
 libraries into 14 representative gst enzymes with known structures 
 and substrates using the plop program [ jacobsonproteins2004 , 55 , 35115048827 ] . 
 for x - ray structures as receptors , 
 one of the top 3 ranked models 
 is within 3  all - atom root mean square deviation ( rmsd ) of the 
 native complex in 11 of the 14 cases ; the enrichment logauc value 
 is better than random in all cases , and better than 25 in 7 of 11 
 cases . for comparative models as receptors , near - native ligand 
  enzyme 
 configurations are often sampled but difficult to rank highly . for 
 models based on templates with the highest sequence identity , 
 the 
 enrichment logauc is better than 25 in 5 of 11 cases , not significantly 
 different from the crystal structures . in conclusion , we show that 
 covalent docking can be a useful tool for substrate discovery and 
 point out specific challenges for future method improvement . 
Introduction
Methods
Results
Discussion
Conclusion