BACKGROUND
molecular docking simulation constitutes one of the main stages of rational or structure-based drug design  <cit> .
it provides a prediction for a molecule binding to a protein in order to form a stable complex  <cit> .
knowledge of proper orientation can be used to predict the strength of association or binding affinity between two molecules.
initially, molecular docking was compared to the classic "key-lock” theory of enzyme-substrate specificity postulated by emil fischer in  <dig> .
in this model, the three-dimensional  structure of both ligand and protein complement each other in the same way a key fits the corresponding lock  <cit> .
however, since both protein and ligand are flexible molecules, the concept is no longer adequate as during the process of molecular docking both ligand and protein adjust their conformation in order to achieve the best protein-ligand fit.
this type of conformational adjustment between the two molecules, or the induced fit theory, was first presented by daniel e. koshland jr. in  <dig>  <cit> .
in order to make molecular docking simulations more realistic, an important issue is to treat both receptor and ligand as flexible structures instead of rigid bodies.
in many methods the ligand, usually a small molecule with up to dozens of atoms, is treated as flexible, but the flexibility of the protein/receptor , depending on their complexity and size, which can reach dozens of thousands of atoms, is still treated in a more restricted manner.
according to cozzini et al.
“the challenge for drug discovery, as well as docking or virtual screening, is to model the plasticity of the receptor so that both structures can conformationally adapt to each other”  <cit> .
therefore, it is well known in the literature that the recognition of the ligand by the receptor is a dynamic event, where both structures change their conformations to minimize the free energy of binding  for their association  <cit> .
nevertheless, most methods of docking employ a single, rigid structure of the receptor.
this happens for practical reasons.
if we try to consider the explicit flexibility of receptor and ligand, the conformational space to be considered quickly becomes impractical  <cit> , as the process would require an enormous computational effort.
in addition, totrov and abagyan  <cit>  state that the best docking algorithms today erroneously predict the position of ligand binding in  <dig> to 70% of the cases, when only one receptor conformation is considered.
in biological systems, proteins express their functions in aqueous or semi-fluid environments.
when in solution, proteins exist in a number of energetically different conformations, so that their structure is best described when all the different states are represented  <cit> .
a set of structures of a particular protein can be determined experimentally by x-ray crystallography or nuclear magnetic resonance, through computational methods which includes monte carlo and molecular dynamics  simulations  <cit> .
therefore if we consider the explicit flexibility using multiple receptor conformations, there are also a number of approaches.
an example is the relaxed complex scheme   <cit> .
the idea is to perform md simulation of the unliganded receptor before docking to address its flexibility.
the rcs method acknowledges that a ligand probably will bind to conformations of the receptor that occur rarely in its dynamic state.
this strong binding often indicates multivalent attachment of the ligand to the receptor.
the second phase of the rcs method involves the rapid docking of small libraries of ligands to a large ensemble of md-derived receptor conformations.
further information and comprehensive reviews of different methods for flexible-receptor flexible-ligand docking can be found in  <cit> .
there have been very few studies on the effect of md-derived receptor flexibility in flexible-ligand docking simulations  <cit> .
in order to improve the body of evidence about the role that receptor flexibility plays in molecular docking simulations we use ensembles of snapshots obtained from three different md simulations of the inha enzyme from mtb; inha_wt, and the mutants inha_i16t, and inha_i21v.
each ensemble of snapshots is denominated a fully-flexible receptor  model.
we have used them to systematically explore their docking to three different inha inhibitors, namely, ethionamide , triclosan , and pentacyanoferrate .
methods
in order to carry out docking simulations, we need a receptor model, and at least one ligand, as well as docking software.
below we describe the different steps involved in processing a docking simulation, including the definition of the rigid and flexible receptor models, the preparation of the ligands, particularly the reference ligands’ position that will be used in our supervised docking simulations.
the single-rigid receptor model
the inha enzyme or 2-trans-enoyl-acp  reductase  from m. tuberculosis was chosen as receptor model for this work because of its importance as a drug target against tuberculosis  <cit> .
it belongs to the sdr  family of proteins, which uses nadh  as coenzyme.
the main feature of this family is the topology of the polypeptide backbone, where each subunit of the protein is composed of a single domain with a classical rossmann-fold topology  <cit> .
it is characterized by a 7-strands parallel β-sheet and eight α-helices, connected by loops and turns, forming the nadh binding site .
the enzyme has a “chair-like” appearance where the “legs” and “backrest” are topologically similar to other dehydrogenases.
the substrate binding cavity is a “pocket” located in the backrest.
it is formed by the substrate-binding loop , strands β <dig>  β <dig>  and helix α <dig>  <cit> .
the nadh coenzyme is positioned in an extended conformation in the pocket along the top of the parallel β-sheet.
the adenine ring is nearly parallel to the “seat” of the structure while the nicotinamide portion faces backwards, pointing to the base of the substrate binding cavity  <cit> .
the a- and b-loops, as well as the substrate-binding loop are structural motifs key to this enzyme function.
the rigid model of inha is taken from its crystal structure  and will be called 1eny throughout the text.
the fully-flexible receptor  models
in this study we have considered three ffr models of mtb’s inha enzyme:  the inha_wt,  the inha_i16t mutant, and  the inha_i21v mutant.
all three ffr models were built from md simulations of the inha_wt-nadh , inha_i16t-nadh  and inha_i21v-nadh  complexes generated with the sander module of amber  <dig>  <cit> .
the mutants were constructed from the  <dig>  ns instantaneous snapshot of the inha_wt md simulation  <cit> .
these mutations have been clinically associated with isoniazid drug resistance  <cit> .
schroeder et al.
<cit>  showed that the resistance was due to a reduced affinity of the nadh for the inha enzyme.
because eth action is similar to that of isoniazid, we expect a difference in its binding to both the inha_wt and inha mutants.
the md trajectories of the mutants and wild type simulations had different lengths.
however, to have coherent flexible models, we considered  <dig>  ns of each md simulation.
for instance, each mutant ffr model was generated from a trajectory in the interval ranging from  <dig>  ns to  <dig>  ns.
we used intervals of  <dig>  ps.
accordingly, each ffr model is composed of  <dig>  snapshots.
receptor preparation
each ffr model of the inha enzyme, which is an ensemble of snapshots from its md simulation, was converted into the conventional pdb file using the ptraj module of amber  <dig>  <cit> .
ptraj also computed an average structure for each ffr model considering the production phase of each trajectory, the last  <dig>  snapshots.
the ffr models were then fitted to their average structures, as well as to the rigid model 1eny.
in this manner, all snapshots making up each ffr model, the average structure and the rigid model of inha will be all in the same frame of reference.
finally, we added the appropriate partial atomic charges, and solvation parameters using the addsol module of autodock  <dig> .
<dig>  <cit> , to each receptor model.
all the steps described above were accomplished with the scientific workflow described in  <cit> .
the ligands
ethionamide   is a relatively small molecule, composed of  <dig> atoms .
this is a powerful second line tuberculostatic, an isoniazid  structural analogue, and is widely used in the treatment of tuberculosis because its primary target is the inha protein.
like inh, eth is also a pro-drug that requires prior activation  <cit> .
its mode of action is similar to that of inh.
eth binds covalently to carbon  <dig> of the nicotinamide portion of nadh to form the adduct eth-nad .
triclosan ,  is composed of  <dig> atoms grouped into two aromatic rings .
it is an antibacterial and antifungal agent commonly found in various preparations ranging from toothpaste, cosmetics in general, antiseptic soaps and even plastic.
in  <dig>  mcmurry, oethinger and levy  <cit>  suggested for the first time that tcl blocked the biosynthesis of fatty acids by inhibiting the enoyl reductase  or inha.
the tcl phenolic ring  forms π-stacking interactions with the nicotinamide ring of nadh.
such interactions are formed due to stacking of aromatic rings of different molecules through van der waals forces  <cit> .
pentacyanoferrate  is the result of a rational drug design effort by santos, basso and co-workers  <cit>  in an attempt to find new inhibitors for m. tuberculosis’ inha enzyme that do not require activation.
this is an inh molecule with a metallic centre, the pentacyanoferrate group, bound to it .
the pif molecule is composed of  <dig> atoms.
since the crystal structure of the inha-pif complex is yet not available, we performed molecular docking simulations to predict the binding mode of pif to inha  <cit> .
definition of the reference ligands
in many docking simulations described in the literature, the approach chosen is blind docking, where the ligand is placed at an initial arbitrary position within the active site of the target receptor, and from there, the docking software seeks to find the best ligand orientation that should correspond to the most negative feb. the docking results are presented in histograms .
the problem is that the most negative estimated feb does not necessarily correspond to the actual binding mode  of the ligand found in experimental results  <cit> .
autodock clusters by first sorting all the docked conformations from lowest  to highest energy .
the best overall docked conformation is used as the ‘seed’ for the first cluster .
then the coordinates of the second best conformation are compared with those of the best to calculate the root-mean square deviation  between the two conformations.
if the calculated rmsd value is smaller than or equal the specified cut-off, which is  <dig>  by default, that conformation is added to the ‘bin’ containing the best conformation.
if not, the second becomes the reference for a second ‘bin’ .
in this example we have four “bins” or clusters.
as can be seen from the table, the  <dig> “runs” generated four clusters.
cluster  <dig> contains six out of  <dig> “runs” and so forth.
the last column shows the rmsd of the corresponding conformation with respect to a reference ligand position.
the bold value is the best estimated docked energy corresponding to the best docked pose expected from experiments  <cit> .
for this reason, before starting the docking simulations with the three different receptor models, we performed a blind docking simulation with the 1eny structure and eth in order to compare the best docking results with the experimental one from the crystal structure  <cit> .
in fact, run  <dig>  is the one that gives the best ligand pose of eth in inha.
note that its feb value of - <dig>  kcal/mol is greater than the best one  ranked automatically by the docking program.
from this test experiment, we decided not to perform blind docking, but supervised dockings instead.
for each ligand, we first tried to find its conformation and orientation corresponding to the best pose.
the coordinates for this ligand pose was then saved and used as a reference coordinate to calculate the reference root-mean-square deviation   of the ligand in the docking simulations.
the same process was repeated for tcl, and pif.
thus, in the docking experiments in this study both, the feb estimates and rmsd values will be important to describe the effect of the ffr models in docking simulations.
ligand preparation
using the vmd software  <cit> , the ligands in the tripos mol <dig> file format were randomly positioned at the receptor binding site in the rigid model, as well as in the three inha ffr models.
since all models are in the same frame of reference, the initial position of all three ligands is same in all docking simulations for all models.
the ligand flexibility was determined by the deftors module of autodock  <dig> .
<dig>  eth and tcl, and pif have two and three rotatable bonds, respectively.
molecular docking simulations
the docking simulations were performed with autodock  <dig> .
<dig>  <cit>  which has been extensively tested and proved to be successful in a variety of docking experiments  <cit> .
it can use various techniques to explore the different conformations a ligand can assume, combining the advantages of a complete search space and the assessment of the best feb  <cit> .
docking simulation parameters
with the rigid 1eny model and the three ffr models of inha in the pdbqs format, the active site was defined within a grid of  <dig> x  <dig> x  <dig> points, spaced at  <dig>  Å. this set up generated a grid box with approximately  <dig> x  <dig> x  <dig> Å <dig>  centred in the initial position of the ligand.
this grid is large enough to include the nadh, as well as the substrate-binding cavity of inha.
each docking simulation was composed of  <dig> independent runs, for which a maximum number of  <dig>  generations were produced, employing the lamarckian genetic algorithm  on the initial population of  <dig> individuals, a maximum number of  <dig>  energy evaluations, with an elitism value of  <dig>  a mutation rate of  <dig> , and a cross-over rate of  <dig> .
for the local search the pseudo-solis and wets method was applied using default parameters.
each run provides one predicted binding mode.
at the end of the docking experiment binding modes with rmsd of  <dig>  Å within each other were placed in the same cluster.
control docking
we performed docking simulations as a control for each of the ligands using the 1eny rigid model.
all results found for the docked ffr models of inha will be compared to these controls.
receptor-ligand interaction analyses
after performing the docking simulations, we employed ligplot  <dig> .
<dig>  <cit>  to analyze the hydrogen bonds and hydrophobic contacts between the ligands and the rigid and ffr models of inha_wt and the mutants inha_i16t and inha_ i21v.
ligplot defines a ffr residue to be in a hydrophobic contact  with a ligand if there is at least one heavy atom of the residue within  <dig>  Å of some atom from the ligand.
for hydrogen bonds  the criteria is more restricted.
the donor and acceptor atoms of ffrs’ residues and ligands have to be at a maximum distance of  <dig>  Å  <cit> .
automating the molecular docking simulations
we carried out the molecular docking experiments in a core  <dig> quad  <dig>  ghz machine, with  <dig> gb of ram and  <dig> gb hd.
however, as there is currently no reliable, automated way to perform docking simulations in ffr models as the ones used in this work, our solution was to create in-house processing scripts using the programming languages bash, awk and python.
automating the docking analyses
in order to carry out the docking analyses with ligplot  <cit> , we also created processing scripts, which are detailed below.
<dig>  for each docked snapshot in a given ffr-ligand complex we extracted and stored the best runs according to the lowest feb  and the lowest rmsd in separate tables.
we note that this rmsd measure is calculated with respect to the reference ligand pose .
this produced nine tables .
<dig>  we then ran ligplot for each ffr-ligand complex, processing the output files in order to extract and store the amino acid contacts of each one into a secondary table.
this step produced  <dig> tables .
<dig>  finally, we counted the amino acid contacts, producing extra  <dig> tables for the ffr-ligand complexes describe above.
the goal here was to identify which residues of the ffr models had interacted most with the ligands.
RESULTS
flexible eth, tcl, and pif ligands were docked to three different ffr models of mtb’s inha: inha_wt, inha_i16t and inha_i21v.
each complex gave us a set of docking results composed of the average and standard deviations  of the best feb  and its corresponding rmsd  with respect to the reference pose for the ligand .
we also obtained similar statistics for the set of feb values matching the lowest rmsd with respect to the reference pose for the ligands .
the latter values are expected to represent the docking results for the best ligand pose .
for each set we also calculated minimum and maximum values of the feb and rmsd in order to compare with our control docking simulation.
based on the average and the standard deviations of the feb values, we do not see much difference between the ffr and the rigid 1eny models.
the only exception is for the pif ligand.
the average best feb for the initial pif position was - <dig>  ±  <dig> , - <dig>  ±  <dig> , and - <dig>  ±  <dig>  kcal/mol while for the reference pose, the average feb was - <dig>  ±  <dig> , - <dig>  ±  <dig>  and - <dig>  ±  <dig>  kcal/mol for inha_wt, inha_i16t and inha_i21v, respectively.
the equivalent feb values for the 1eny rigid model was - <dig>  and - <dig>  kcal/mol, respectively.
these values represent a difference that varies from  <dig>  to  <dig>  kcal/mol between the rigid and ffr models of inha.
also, they are greater than the intrinsic error  attributed to the estimation of the feb by autodock <dig> .
<dig>  <cit> .
the explicit flexibility of the inha receptor clearly had an impact in the way pif interacted with it.
the first column indicates the ligands.
the second column describes the statistics evaluated for each docking result: the average, standard deviation , the minimum and maximum feb values.
column  <dig>  contains the values of best feb in kcal/mol and column  <dig>  its related rmsd in Å. columns  <dig>  and  <dig>  show the feb corresponding to the lowest rmsd calculated with respect to the reference ligand position .
columns 3- <dig> contain the results for inha_wt.
columns 7- <dig> and 11- <dig> are the values of a, b, c and d for inha_i16t and inha_i21v mutants, respectively.
the last row shows the results for the 1eny rigid model.
analyses of ffr models-ligands interactions
the coordinates of the ffr models-ligands complexes were analysed by ligplot.
through this analysis we were able to identify which residues were making hhb and nnb contacts with the ligands in at least one of the ffrs models’ snapshots  which means that, in order to identify those residues we simply considered all snapshots making at least one contact, either hhb or nnb, with each ligand.
when the same residue made both intermolecular contacts, the redundancy was eliminated.
as table  <dig> shows, there is a significant increase in the number of residues that are able to interact with the ligands in the ffr models as compared to the rigid one.
overall, a minimum of nine and a maximum of  <dig> different residues interacted with any of the three ligands in all ffr models.
conversely, in the rigid 1eny model these numbers vary from only zero to five.
these data further corroborate our hypothesis that considering the explicit flexibility of the receptor during docking simulations allows a ligand to interact, even casually, with other receptor residues than those found in a single, rigid crystal structure.
the behaviour of the ffr models simulates the induced fit phenomenon  <cit> .
the total number of receptor residues that made hhb and nnb contacts with the ligands and the union of these residues  for the best feb  and for the reference or best pose .
from the data in table  <dig> we were able to single out the top  <dig> ffr amino acid residues that made contact with eth and pif ligands, and the top  <dig> ffr amino acid residues that made contact with tcl, in at least one of their snapshots .
the first columns show, from the total of amino acid residues found in table  <dig>  the top  <dig>  residues that were able to make contacts, hhb and/or nnb, with the ligands  eth,  tcl, and  pif, respectively, in at least 10% of the trajectories.
note that each ffr model is a md simulation trajectory.
columns  <dig> to  <dig> show the number of snapshots for each ffr model involved in such contacts .
note that we are not trying to count the absolute number of contacts.
instead we are counting how many snapshots of each ffr model provided those contacts.
for eth  we found three polar , one acidic and one basic , and  <dig> hydrophobic residues.
as for tcl  we found three polar , one basic , and  <dig> hydrophobic residues.
for pif  we found two polar , one acidic and one basic , and  <dig> apolar residues.
as can be seen in column  <dig> of table  <dig>  most of the residues are hydrophobic, which is somewhat expected since the inha active site is predominantly apolar.
figure  <dig> shows where these residues are located in the receptor structure.
these residues are part of both, the nadh and, in particular, the substrate binding pockets.
discussion
in order to make molecular docking simulations more realistic, an important issue is to treat both receptor and ligand as flexible structures instead of rigid bodies.
however, this has been demonstrated not to be a trivial task  <cit> .
we know that the enzyme inha from m. tuberculosis is a highly flexible receptor  <cit>  and is a very important molecular target for the development of new drugs against tuberculosis  <cit> .
given that only a few studies  <cit>  have addressed the role of the explicit flexibility of receptors in molecular docking simulations and the biomedical importance of the inha enzyme, here we reported our systematic analysis of the effect of inha explicit flexibility in docking simulations to three ligands known to be its inhibitors.
our three ffr models of inha, namely inha_wt, inha_i16t, and inha_i21v were generated from their md simulations trajectories  <cit> .
as a result we obtained a number of 3-d snapshots of the protein that were slightly different from each other.
after that, we then docked each ligand  to each one of ffr models of inha.
we believe that the structural differences from one snapshot to another, creates a space in the receptors binding cavities, which differs from the one we see in the rigid, crystal structure.
for instance, in the case of the inha_wt ffr model, and using the castp server  <cit> , the snapshots at  <dig>  ns and  <dig>  ns have both  <dig>  Å <dig> and  <dig>  Å <dig>  respectively.
there is a variation of approximately  <dig> Å <dig> in the volume of the inha’s major binding pocket between these two snapshots.
this “new space” gives the ligand an opportunity to better explore the binding cavities, increasing its chance to accommodate to it.
due to the fact that we knew beforehand how eth and tcl inhibits inha, blind docking simulations would not be appropriate.
eth binds covalently to carbon  <dig> of the nicotinamide portion of nadh to form the eth-nad adduct  <cit> , whereas the phenolic ring of tcl forms π-stacking interaction with the nadh nicotinamide ring  <cit> .
to acknowledge that, in our docking simulations, although the ligands were randomly positioned in the receptor active site, we created what we called the reference ligand.
this means that the rmsd calculated during the docking simulations did not use the initial position of the ligand, but rather the position of the reference ligand or best pose, closer to the one expected to form the adduct  <cit> .
our analyses showed that up to  <dig> different receptors’ residues interact with the ligands as opposed to up to only five residues in the rigid 1eny model.
this constitutes, in our opinion, a strong basis to recommend the use of explicitly flexible models of mtb’s inha enzyme in virtual screening efforts to search for novel drug candidates against tuberculosis.
CONCLUSIONS
with our data analyses we were able to find a total of up to  <dig> receptor amino acid residues interacting with the ligands employed in this study.
performing docking under the same conditions, but in the rigid, crystal structure 1eny, we were able to find only five for eth and two for both tcl and pif.
these numbers supports our hypothesis that flexible receptor models can accommodate a more diverse range of ligand conformations.
this indicates that they are more prone to select a new ligand capable of binding to inha than they would do if we used only one receptor conformation.
in other words, taking the receptor plasticity into account when performing docking simulation means that amino acid residues, loops and turns, can move slightly in different directions, giving the ligand a better chance to accommodate itself in the receptors’ binding site.
nonetheless, we are aware that our results were based on a short md simulation; only  <dig>  ns long and that nowadays much longer simulations can be generated.
hence, one of our future goals is to consider whether or not longer md simulations would affect our conclusions, or even how our docking simulations would behave if they were performed on different ffr models generated from md simulations at physiological temperatures.
competing interests
the authors declare that they have no competing interests.
authors' contributions
eml, ksm, and ons conceived and designed the experiments.
eml performed the experiments, analyzed the data and wrote the manuscript first draft.
mc wrote and executed the automating processing scripts .
ons wrote the final version of the manuscript.
all authors read and approved the final version of the manuscript.

 molecular docking simulation constitutes one of the main stages of rational or structure-based drug design. it provides a prediction for a molecule binding to a protein in order to form a stable complex. knowledge of proper orientation can be used to predict the strength of association or binding affinity between two molecules. initially, molecular docking was compared to the classic "key-lock” theory of enzyme-substrate specificity postulated by emil fischer in. in this model, the three-dimensional structure of both ligand and protein complement each other in the same way a key fits the corresponding lock. however, since both protein and ligand are flexible molecules, the concept is no longer adequate as during the process of molecular docking both ligand and protein adjust their conformation in order to achieve the best protein-ligand fit. this type of conformational adjustment between the two molecules, or the induced fit theory, was first presented by daniel e. koshland jr. in. in order to make molecular docking simulations more realistic, an important issue is to treat both receptor and ligand as flexible structures instead of rigid bodies. in many methods the ligand, usually a small molecule with up to dozens of atoms, is treated as flexible, but the flexibility of the protein/receptor, depending on their complexity and size, which can reach dozens of thousands of atoms, is still treated in a more restricted manner. according to cozzini et al. “the challenge for drug discovery, as well as docking or virtual screening, is to model the plasticity of the receptor so that both structures can conformationally adapt to each other”. therefore, it is well known in the literature that the recognition of the ligand by the receptor is a dynamic event, where both structures change their conformations to minimize the free energy of binding for their association. nevertheless, most methods of docking employ a single, rigid structure of the receptor. this happens for practical reasons. if we try to consider the explicit flexibility of receptor and ligand, the conformational space to be considered quickly becomes impractical, as the process would require an enormous computational effort. in addition, totrov and abagyan state that the best docking algorithms today erroneously predict the position of ligand binding in to 70% of the cases, when only one receptor conformation is considered. in biological systems, proteins express their functions in aqueous or semi-fluid environments. when in solution, proteins exist in a number of energetically different conformations, so that their structure is best described when all the different states are represented. a set of structures of a particular protein can be determined experimentally by x-ray crystallography or nuclear magnetic resonance, through computational methods which includes monte carlo and molecular dynamics simulations. therefore if we consider the explicit flexibility using multiple receptor conformations, there are also a number of approaches. an example is the relaxed complex scheme. the idea is to perform md simulation of the unliganded receptor before docking to address its flexibility. the rcs method acknowledges that a ligand probably will bind to conformations of the receptor that occur rarely in its dynamic state. this strong binding often indicates multivalent attachment of the ligand to the receptor. the second phase of the rcs method involves the rapid docking of small libraries of ligands to a large ensemble of md-derived receptor conformations. further information and comprehensive reviews of different methods for flexible-receptor flexible-ligand docking can be found in. there have been very few studies on the effect of md-derived receptor flexibility in flexible-ligand docking simulations. in order to improve the body of evidence about the role that receptor flexibility plays in molecular docking simulations we use ensembles of snapshots obtained from three different md simulations of the inha enzyme from mtb; inha_wt, and the mutants inha_i16t, and inha_i21v. each ensemble of snapshots is denominated a fully-flexible receptor model. we have used them to systematically explore their docking to three different inha inhibitors, namely, ethionamide, triclosan, and pentacyanoferrate. methods in order to carry out docking simulations, we need a receptor model, and at least one ligand, as well as docking software. below we describe the different steps involved in processing a docking simulation, including the definition of the rigid and flexible receptor models, the preparation of the ligands, particularly the reference ligands’ position that will be used in our supervised docking simulations. the single-rigid receptor model the inha enzyme or 2-trans-enoyl-acp reductase from m. tuberculosis was chosen as receptor model for this work because of its importance as a drug target against tuberculosis. it belongs to the sdr family of proteins, which uses nadh as coenzyme. the main feature of this family is the topology of the polypeptide backbone, where each subunit of the protein is composed of a single domain with a classical rossmann-fold topology. it is characterized by a 7-strands parallel β-sheet and eight α-helices, connected by loops and turns, forming the nadh binding site. the enzyme has a “chair-like” appearance where the “legs” and “backrest” are topologically similar to other dehydrogenases. the substrate binding cavity is a “pocket” located in the backrest. it is formed by the substrate-binding loop, strands β β and helix α. the nadh coenzyme is positioned in an extended conformation in the pocket along the top of the parallel β-sheet. the adenine ring is nearly parallel to the “seat” of the structure while the nicotinamide portion faces backwards, pointing to the base of the substrate binding cavity. the a- and b-loops, as well as the substrate-binding loop are structural motifs key to this enzyme function. the rigid model of inha is taken from its crystal structure and will be called 1eny throughout the text. the fully-flexible receptor models in this study we have considered three ffr models of mtb’s inha enzyme: the inha_wt, the inha_i16t mutant, and the inha_i21v mutant. all three ffr models were built from md simulations of the inha_wt-nadh, inha_i16t-nadh and inha_i21v-nadh complexes generated with the sander module of amber. the mutants were constructed from the ns instantaneous snapshot of the inha_wt md simulation. these mutations have been clinically associated with isoniazid drug resistance. schroeder et al. showed that the resistance was due to a reduced affinity of the nadh for the inha enzyme. because eth action is similar to that of isoniazid, we expect a difference in its binding to both the inha_wt and inha mutants. the md trajectories of the mutants and wild type simulations had different lengths. however, to have coherent flexible models, we considered ns of each md simulation. for instance, each mutant ffr model was generated from a trajectory in the interval ranging from ns to ns. we used intervals of ps. accordingly, each ffr model is composed of snapshots. receptor preparation each ffr model of the inha enzyme, which is an ensemble of snapshots from its md simulation, was converted into the conventional pdb file using the ptraj module of amber. ptraj also computed an average structure for each ffr model considering the production phase of each trajectory, the last snapshots. the ffr models were then fitted to their average structures, as well as to the rigid model 1eny. in this manner, all snapshots making up each ffr model, the average structure and the rigid model of inha will be all in the same frame of reference. finally, we added the appropriate partial atomic charges, and solvation parameters using the addsol module of autodock., to each receptor model. all the steps described above were accomplished with the scientific workflow described in. the ligands ethionamide is a relatively small molecule, composed of atoms. this is a powerful second line tuberculostatic, an isoniazid structural analogue, and is widely used in the treatment of tuberculosis because its primary target is the inha protein. like inh, eth is also a pro-drug that requires prior activation. its mode of action is similar to that of inh. eth binds covalently to carbon of the nicotinamide portion of nadh to form the adduct eth-nad. triclosan, is composed of atoms grouped into two aromatic rings. it is an antibacterial and antifungal agent commonly found in various preparations ranging from toothpaste, cosmetics in general, antiseptic soaps and even plastic. in mcmurry, oethinger and levy suggested for the first time that tcl blocked the biosynthesis of fatty acids by inhibiting the enoyl reductase or inha. the tcl phenolic ring forms π-stacking interactions with the nicotinamide ring of nadh. such interactions are formed due to stacking of aromatic rings of different molecules through van der waals forces. pentacyanoferrate is the result of a rational drug design effort by santos, basso and co-workers in an attempt to find new inhibitors for m. tuberculosis’ inha enzyme that do not require activation. this is an inh molecule with a metallic centre, the pentacyanoferrate group, bound to it. the pif molecule is composed of atoms. since the crystal structure of the inha-pif complex is yet not available, we performed molecular docking simulations to predict the binding mode of pif to inha. definition of the reference ligands in many docking simulations described in the literature, the approach chosen is blind docking, where the ligand is placed at an initial arbitrary position within the active site of the target receptor, and from there, the docking software seeks to find the best ligand orientation that should correspond to the most negative feb. the docking results are presented in histograms. the problem is that the most negative estimated feb does not necessarily correspond to the actual binding mode of the ligand found in experimental results. autodock clusters by first sorting all the docked conformations from lowest to highest energy. the best overall docked conformation is used as the ‘seed’ for the first cluster. then the coordinates of the second best conformation are compared with those of the best to calculate the root-mean square deviation between the two conformations. if the calculated rmsd value is smaller than or equal the specified cut-off, which is by default, that conformation is added to the ‘bin’ containing the best conformation. if not, the second becomes the reference for a second ‘bin’. in this example we have four “bins” or clusters. as can be seen from the table, the “runs” generated four clusters. cluster contains six out of “runs” and so forth. the last column shows the rmsd of the corresponding conformation with respect to a reference ligand position. the bold value is the best estimated docked energy corresponding to the best docked pose expected from experiments. for this reason, before starting the docking simulations with the three different receptor models, we performed a blind docking simulation with the 1eny structure and eth in order to compare the best docking results with the experimental one from the crystal structure. in fact, run is the one that gives the best ligand pose of eth in inha. note that its feb value of - kcal/mol is greater than the best one ranked automatically by the docking program. from this test experiment, we decided not to perform blind docking, but supervised dockings instead. for each ligand, we first tried to find its conformation and orientation corresponding to the best pose. the coordinates for this ligand pose was then saved and used as a reference coordinate to calculate the reference root-mean-square deviation of the ligand in the docking simulations. the same process was repeated for tcl, and pif. thus, in the docking experiments in this study both, the feb estimates and rmsd values will be important to describe the effect of the ffr models in docking simulations. ligand preparation using the vmd software, the ligands in the tripos mol file format were randomly positioned at the receptor binding site in the rigid model, as well as in the three inha ffr models. since all models are in the same frame of reference, the initial position of all three ligands is same in all docking simulations for all models. the ligand flexibility was determined by the deftors module of autodock. eth and tcl, and pif have two and three rotatable bonds, respectively. molecular docking simulations the docking simulations were performed with autodock. which has been extensively tested and proved to be successful in a variety of docking experiments. it can use various techniques to explore the different conformations a ligand can assume, combining the advantages of a complete search space and the assessment of the best feb. docking simulation parameters with the rigid 1eny model and the three ffr models of inha in the pdbqs format, the active site was defined within a grid of x x points, spaced at Å. this set up generated a grid box with approximately x x Å centred in the initial position of the ligand. this grid is large enough to include the nadh, as well as the substrate-binding cavity of inha. each docking simulation was composed of independent runs, for which a maximum number of generations were produced, employing the lamarckian genetic algorithm on the initial population of individuals, a maximum number of energy evaluations, with an elitism value of a mutation rate of, and a cross-over rate of. for the local search the pseudo-solis and wets method was applied using default parameters. each run provides one predicted binding mode. at the end of the docking experiment binding modes with rmsd of Å within each other were placed in the same cluster. control docking we performed docking simulations as a control for each of the ligands using the 1eny rigid model. all results found for the docked ffr models of inha will be compared to these controls. receptor-ligand interaction analyses after performing the docking simulations, we employed ligplot. to analyze the hydrogen bonds and hydrophobic contacts between the ligands and the rigid and ffr models of inha_wt and the mutants inha_i16t and inha_ i21v. ligplot defines a ffr residue to be in a hydrophobic contact with a ligand if there is at least one heavy atom of the residue within Å of some atom from the ligand. for hydrogen bonds the criteria is more restricted. the donor and acceptor atoms of ffrs’ residues and ligands have to be at a maximum distance of Å. automating the molecular docking simulations we carried out the molecular docking experiments in a core quad ghz machine, with gb of ram and gb hd. however, as there is currently no reliable, automated way to perform docking simulations in ffr models as the ones used in this work, our solution was to create in-house processing scripts using the programming languages bash, awk and python. automating the docking analyses in order to carry out the docking analyses with ligplot, we also created processing scripts, which are detailed below. for each docked snapshot in a given ffr-ligand complex we extracted and stored the best runs according to the lowest feb and the lowest rmsd in separate tables. we note that this rmsd measure is calculated with respect to the reference ligand pose. this produced nine tables. we then ran ligplot for each ffr-ligand complex, processing the output files in order to extract and store the amino acid contacts of each one into a secondary table. this step produced tables. finally, we counted the amino acid contacts, producing extra tables for the ffr-ligand complexes describe above. the goal here was to identify which residues of the ffr models had interacted most with the ligands.  flexible eth, tcl, and pif ligands were docked to three different ffr models of mtb’s inha: inha_wt, inha_i16t and inha_i21v. each complex gave us a set of docking results composed of the average and standard deviations of the best feb and its corresponding rmsd with respect to the reference pose for the ligand. we also obtained similar statistics for the set of feb values matching the lowest rmsd with respect to the reference pose for the ligands. the latter values are expected to represent the docking results for the best ligand pose. for each set we also calculated minimum and maximum values of the feb and rmsd in order to compare with our control docking simulation. based on the average and the standard deviations of the feb values, we do not see much difference between the ffr and the rigid 1eny models. the only exception is for the pif ligand. the average best feb for the initial pif position was - ±, - ±, and - ± kcal/mol while for the reference pose, the average feb was - ±, - ± and - ± kcal/mol for inha_wt, inha_i16t and inha_i21v, respectively. the equivalent feb values for the 1eny rigid model was - and - kcal/mol, respectively. these values represent a difference that varies from to kcal/mol between the rigid and ffr models of inha. also, they are greater than the intrinsic error attributed to the estimation of the feb by autodock.. the explicit flexibility of the inha receptor clearly had an impact in the way pif interacted with it. the first column indicates the ligands. the second column describes the statistics evaluated for each docking result: the average, standard deviation, the minimum and maximum feb values. column contains the values of best feb in kcal/mol and column its related rmsd in Å. columns and show the feb corresponding to the lowest rmsd calculated with respect to the reference ligand position. columns 3- contain the results for inha_wt. columns 7- and 11- are the values of a, b, c and d for inha_i16t and inha_i21v mutants, respectively. the last row shows the results for the 1eny rigid model. analyses of ffr models-ligands interactions the coordinates of the ffr models-ligands complexes were analysed by ligplot. through this analysis we were able to identify which residues were making hhb and nnb contacts with the ligands in at least one of the ffrs models’ snapshots which means that, in order to identify those residues we simply considered all snapshots making at least one contact, either hhb or nnb, with each ligand. when the same residue made both intermolecular contacts, the redundancy was eliminated. as table shows, there is a significant increase in the number of residues that are able to interact with the ligands in the ffr models as compared to the rigid one. overall, a minimum of nine and a maximum of different residues interacted with any of the three ligands in all ffr models. conversely, in the rigid 1eny model these numbers vary from only zero to five. these data further corroborate our hypothesis that considering the explicit flexibility of the receptor during docking simulations allows a ligand to interact, even casually, with other receptor residues than those found in a single, rigid crystal structure. the behaviour of the ffr models simulates the induced fit phenomenon. the total number of receptor residues that made hhb and nnb contacts with the ligands and the union of these residues for the best feb and for the reference or best pose. from the data in table we were able to single out the top ffr amino acid residues that made contact with eth and pif ligands, and the top ffr amino acid residues that made contact with tcl, in at least one of their snapshots. the first columns show, from the total of amino acid residues found in table the top residues that were able to make contacts, hhb and/or nnb, with the ligands eth, tcl, and pif, respectively, in at least 10% of the trajectories. note that each ffr model is a md simulation trajectory. columns to show the number of snapshots for each ffr model involved in such contacts. note that we are not trying to count the absolute number of contacts. instead we are counting how many snapshots of each ffr model provided those contacts. for eth we found three polar, one acidic and one basic, and hydrophobic residues. as for tcl we found three polar, one basic, and hydrophobic residues. for pif we found two polar, one acidic and one basic, and apolar residues. as can be seen in column of table most of the residues are hydrophobic, which is somewhat expected since the inha active site is predominantly apolar. figure shows where these residues are located in the receptor structure. these residues are part of both, the nadh and, in particular, the substrate binding pockets. discussion in order to make molecular docking simulations more realistic, an important issue is to treat both receptor and ligand as flexible structures instead of rigid bodies. however, this has been demonstrated not to be a trivial task. we know that the enzyme inha from m. tuberculosis is a highly flexible receptor and is a very important molecular target for the development of new drugs against tuberculosis. given that only a few studies have addressed the role of the explicit flexibility of receptors in molecular docking simulations and the biomedical importance of the inha enzyme, here we reported our systematic analysis of the effect of inha explicit flexibility in docking simulations to three ligands known to be its inhibitors. our three ffr models of inha, namely inha_wt, inha_i16t, and inha_i21v were generated from their md simulations trajectories. as a result we obtained a number of 3-d snapshots of the protein that were slightly different from each other. after that, we then docked each ligand to each one of ffr models of inha. we believe that the structural differences from one snapshot to another, creates a space in the receptors binding cavities, which differs from the one we see in the rigid, crystal structure. for instance, in the case of the inha_wt ffr model, and using the castp server, the snapshots at ns and ns have both Å and Å respectively. there is a variation of approximately Å in the volume of the inha’s major binding pocket between these two snapshots. this “new space” gives the ligand an opportunity to better explore the binding cavities, increasing its chance to accommodate to it. due to the fact that we knew beforehand how eth and tcl inhibits inha, blind docking simulations would not be appropriate. eth binds covalently to carbon of the nicotinamide portion of nadh to form the eth-nad adduct, whereas the phenolic ring of tcl forms π-stacking interaction with the nadh nicotinamide ring. to acknowledge that, in our docking simulations, although the ligands were randomly positioned in the receptor active site, we created what we called the reference ligand. this means that the rmsd calculated during the docking simulations did not use the initial position of the ligand, but rather the position of the reference ligand or best pose, closer to the one expected to form the adduct. our analyses showed that up to different receptors’ residues interact with the ligands as opposed to up to only five residues in the rigid 1eny model. this constitutes, in our opinion, a strong basis to recommend the use of explicitly flexible models of mtb’s inha enzyme in virtual screening efforts to search for novel drug candidates against tuberculosis.  with our data analyses we were able to find a total of up to receptor amino acid residues interacting with the ligands employed in this study. performing docking under the same conditions, but in the rigid, crystal structure 1eny, we were able to find only five for eth and two for both tcl and pif. these numbers supports our hypothesis that flexible receptor models can accommodate a more diverse range of ligand conformations. this indicates that they are more prone to select a new ligand capable of binding to inha than they would do if we used only one receptor conformation. in other words, taking the receptor plasticity into account when performing docking simulation means that amino acid residues, loops and turns, can move slightly in different directions, giving the ligand a better chance to accommodate itself in the receptors’ binding site. nonetheless, we are aware that our results were based on a short md simulation; only ns long and that nowadays much longer simulations can be generated. hence, one of our future goals is to consider whether or not longer md simulations would affect our conclusions, or even how our docking simulations would behave if they were performed on different ffr models generated from md simulations at physiological temperatures. competing interests the authors declare that they have no competing interests. authors' contributions eml, ksm, and ons conceived and designed the experiments. eml performed the experiments, analyzed the data and wrote the manuscript first draft. mc wrote and executed the automating processing scripts. ons wrote the final version of the manuscript. all authors read and approved the final version of the manuscript. 
protein/receptor explicit flexibility has recently become an important feature of molecular docking simulations. we hereby use ensembles of snapshots obtained from three different md simulations of the inha enzyme from m. tuberculosis, the wild-type, inha_i16t, and inha_i21v mutants to model their explicit flexibility, and to systematically explore their effect in docking simulations with three different inha inhibitors, namely, ethionamide, triclosan, and pentacyanoferrate. the comparison between the rigid crystal structure and the ffr models showed that the inclusion of explicit flexibility, despite the limitations of the ffr models employed in this study, accounts in a substantial manner to the induced fit expected when a protein/receptor and ligand approach each other to interact in the most favourable manner. protein/receptor explicit flexibility, or ffr models, represented as an ensemble of md simulation snapshots, can lead to a more realistic representation of the induced fit effect expected in the encounter and proper docking of receptors to ligands. conversely, in the ffr models this number grows up to an astonishing different residues. the ffr models of inha explicitly characterizes the overall movements of the amino acid residues in helices, strands, loops, and turns, allowing the ligand to properly accommodate itself in the receptor’s binding site. we have produced a short video sequence of each ligand docked to the ffr models of inha_wt. the use of fully-flexible receptor models of inha_wt, inha_i16t, and inha_i21v mutants in docking simulation with the inhibitors eth, tcl, and pif revealed significant differences in the way they interact as compared to the rigid, inha crystal structure. utilization of the intrinsic flexibility of mtb’s inha enzyme and its mutants in virtual screening via molecular docking simulation may provide a novel platform to guide the rational or dynamical-structure-based drug design of novel inhibitors for mtb’s inha.
 protein/receptor explicit flexibility has recently become an important feature of molecular docking simulations. taking the flexibility into account brings the docking simulation closer to the receptors’ real behaviour in its natural environment. several approaches have been developed to address this problem. among them, modelling the full flexibility as an ensemble of snapshots derived from a molecular dynamics simulation of the receptor has proved very promising. despite its potential, however, only a few studies have employed this method to probe its effect in molecular docking simulations. we hereby use ensembles of snapshots obtained from three different md simulations of the inha enzyme from m. tuberculosis, the wild-type, inha_i16t, and inha_i21v mutants to model their explicit flexibility, and to systematically explore their effect in docking simulations with three different inha inhibitors, namely, ethionamide, triclosan, and pentacyanoferrate.  the use of fully-flexible receptor models of inha_wt, inha_i16t, and inha_i21v mutants in docking simulation with the inhibitors eth, tcl, and pif revealed significant differences in the way they interact as compared to the rigid, inha crystal structure. in the latter, only up to five receptor residues interact with the three different ligands. conversely, in the ffr models this number grows up to an astonishing different residues. the comparison between the rigid crystal structure and the ffr models showed that the inclusion of explicit flexibility, despite the limitations of the ffr models employed in this study, accounts in a substantial manner to the induced fit expected when a protein/receptor and ligand approach each other to interact in the most favourable manner.  protein/receptor explicit flexibility, or ffr models, represented as an ensemble of md simulation snapshots, can lead to a more realistic representation of the induced fit effect expected in the encounter and proper docking of receptors to ligands. the ffr models of inha explicitly characterizes the overall movements of the amino acid residues in helices, strands, loops, and turns, allowing the ligand to properly accommodate itself in the receptor’s binding site. utilization of the intrinsic flexibility of mtb’s inha enzyme and its mutants in virtual screening via molecular docking simulation may provide a novel platform to guide the rational or dynamical-structure-based drug design of novel inhibitors for mtb’s inha. we have produced a short video sequence of each ligand docked to the ffr models of inha_wt. these videos are available at http://www.inf.pucrs.br/~osmarns/labio/videos_cohen_et_al_19_07_ htm. 15- november 6th international conference of the brazilian association for bioinformatics and computational biology ouro preto, brazil 