Email updates

Keep up to date with the latest news and content from BMC Genomics and BioMed Central.

This article is part of the supplement: Proceedings of the 6th International Conference of the Brazilian Association for Bioinformatics and Computational Biology (X-meeting 2010)

Open Access Proceedings

Effect of the explicit flexibility of the InhA enzyme from Mycobacterium tuberculosis in molecular docking simulations

Elisangela ML Cohen12, Karina S Machado1, Marcelo Cohen1 and Osmar Norberto de Souza12*

Author Affiliations

1 Laboratório de Bioinformática, Modelagem e Simulação de Biossistemas - LABIO, Programa de Pós-Graduação em Ciência da Computação - PPGCC - Faculdade de Informática - PUCRS - Av. Ipiranga, 6681 Prédio 32 - Sala 608. CEP: 90619-900, Porto Alegre, RS, Brasil

2 Programa de Pós-Graduação em Biologia Celular e Molecular - PPGBCM - Faculdade de Biociências - PUCRS - Av. Ipiranga, 6681 Prédio 12. CEP: 90619-900, Porto Alegre, RS, Brasil

For all author emails, please log on.

BMC Genomics 2011, 12(Suppl 4):S7  doi:10.1186/1471-2164-12-S4-S7

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2164/12/S4/S7


Published:22 December 2011

© 2011 Cohen et al; licensee BioMed Central Ltd.

This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

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 (MD) 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 (Mtb), the wild-type (InhA_wt), 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 (ETH), triclosan (TCL), and pentacyano(isoniazid)ferrate(II) (PIF).

Results

The use of fully-flexible receptor (FFR) 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 (PDB ID: 1ENY). 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 80 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.

Conclusions

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 (ETH, TCL and PIF) 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_2011.htm webcite.

Background

Molecular docking simulation constitutes one of the main stages of rational or structure-based drug design [1]. It provides a prediction for a molecule binding to a protein in order to form a stable complex [2]. 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 1894 (Reviewed by Koshland Jr., [3,4]). In this model, the three-dimensional (3-D) structure of both ligand and protein complement each other in the same way a key fits the corresponding lock [5]. 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 1958 [3,4].

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 (for simplicity, herein protein and receptor are synonymous), 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” [6]. 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 (FEB) for their association [7]. 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 [8,9], as the process would require an enormous computational effort. In addition, Totrov and Abagyan [10] state that the best docking algorithms today erroneously predict the position of ligand binding in 50 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 [6]. 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 (MD) simulations [11]. 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 (RCS) [12]. 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 [13-17].

There have been very few studies on the effect of MD-derived receptor flexibility in flexible-ligand docking simulations [14,16-18]. 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 (FFR) model. We have used them to systematically explore their docking to three different InhA inhibitors, namely, ethionamide (ETH), triclosan (TCL), and pentacyano(isoniazid)ferrate(II) (PIF).

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 (CoA) reductase (EC number 1.3.1.9) from M. tuberculosis was chosen as receptor model for this work because of its importance as a drug target against tuberculosis [19]. It belongs to the SDR (short chain dehydrogenase / reductase) family of proteins, which uses NADH (β-nicotinamide adenine dinucleotide, reduced form) 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 [20]. It is characterized by a 7-strands parallel β-sheet and eight α-helices, connected by loops and turns, forming the NADH binding site (Figure 1). 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 (helices α6 and α7 in Figure 1), strands β4, β5, and helix α5 [20]. 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 [20]. 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 (PDB ID: 1ENY) and will be called 1ENY throughout the text.

thumbnailFigure 1. The InhA tertiary structure. Ribbons representation of Mtb’s InhA (PDB ID:1ENY) rigid receptor model (colored by secondary structure) in complex with the coenzyme NADH (in metallic grey). In yellow are the 7-strands parallel β-sheet and in magenta the eight α-helices, connected by loops (in cyan) and turns (in white). The protein’s N-terminus is composed by helices α1 and α2, and by the β1 to β3 strands, while the C-terminus is formed by helices α7 and α8. Figure produced with VMD [29].

The Fully-Flexible Receptor (FFR) models

In this study we have considered three FFR models of Mtb’s InhA enzyme: (1) the InhA_wt, (2) the InhA_I16T mutant, and (3) the InhA_I21V mutant. All three FFR models were built from MD simulations of the InhA_wt-NADH (3.1 ns), InhA_I16T-NADH (5.5 ns) and InhA_I21V-NADH (6.5 ns) complexes generated with the SANDER module of AMBER 6 [21]. The mutants were constructed from the 1.9 ns instantaneous snapshot of the InhA_wt MD simulation [21]. These mutations have been clinically associated with isoniazid drug resistance [22]. Schroeder et al. [21] 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 3.1 ns of each MD simulation. For instance, each mutant FFR model was generated from a trajectory in the interval ranging from 1.9 ns to 5.0 ns. We used intervals of 1.0 ps. Accordingly, each FFR model is composed of 3,100 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 9 [23]. Ptraj also computed an average structure for each FFR model considering the production phase of each trajectory, the last 1,000 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 3.0.5 [24], to each receptor model. All the steps described above were accomplished with the scientific workflow described in [25].

The ligands

Ethionamide (ETH) (ZINC code: 3872520, accessed on 08/10/2010) is a relatively small molecule, composed of 21 atoms (Figure 2a). This is a powerful second line tuberculostatic, an isoniazid (INH) 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 [26]. Its mode of action is similar to that of INH. ETH binds covalently to carbon 4 of the nicotinamide portion of NADH to form the adduct ETH-NAD (Figure 2a).

thumbnailFigure 2. The ligands. (A) On the left, stick model representation of the ETH ligand. The hydrogen, nitrogen, carbon, and sulphur atoms are colored white, dark blue, cyan, and light brown, respectively. In the right, stick model representation of the adduct ETH-NAD (PDB ID: 2H9I). Activated ETH is colored in yellow and NAD in metallic grey. The adduct binds to the InhA active site with ETH blocking part of the substrate binding cavity. (B) Stick model representation of the TCL ligand. The color scheme and receptor representation are the same as in (A), except for oxygen and chlorine which are colored red and magenta, respectively. Stick model representation of the π stacking interaction between TCL’s A ring (dark blue) and the nicotinamide ring of NADH (from PDB ID: 1P45). (C) Stick model representation of the PIF ligand. Here the iron atom is colored green. Stick model representation of one possible interaction between PIF (green) and the InhA enzyme. InhA main-chain is represented by ribbons.

Triclosan (TCL), (ZINC code: 2216, accessed on 08/10/2010) is composed of 24 atoms grouped into two aromatic rings (Figure 2b). It is an antibacterial and antifungal agent commonly found in various preparations ranging from toothpaste, cosmetics in general, antiseptic soaps and even plastic. In 1998, McMurry, Oethinger and Levy [27] suggested for the first time that TCL blocked the biosynthesis of fatty acids by inhibiting the enoyl reductase (ENR) or InhA. The TCL phenolic ring (A ring in Figure 2b) 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 [28].

Pentacyano(isoniazid)ferrate(II) (PIF) is the result of a rational drug design effort by Santos, Basso and co-workers [29] 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 (Figure 2c). The PIF molecule is composed of 28 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 [29].

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 (Table 1). The problem is that the most negative estimated FEB does not necessarily correspond to the actual binding mode (here we call Best Pose) of the ligand found in experimental results [26].

Table 1. A typical example of a RMSD (in Å) clustering analysis of a 10 “runs” docking simulation with AutoDock 3.0.5

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 [26]. In fact, run 7 (Table 1, in bold) is the one that gives the best ligand pose of ETH in InhA. Note that its FEB value of -8.27 kcal/mol is greater than the best one (-8.99 kcal/mol) 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 (RMSD) (Table 1) 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 [31], the ligands in the Tripos Mol2 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 3.0.5. ETH and TCL, and PIF have two and three rotatable bonds, respectively.

Molecular docking simulations

The docking simulations were performed with AutoDock 3.0.5 [32] which has been extensively tested and proved to be successful in a variety of docking experiments [14]. 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 [32].

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 100 x 60 x 60 points, spaced at 0.375 Å. This set up generated a grid box with approximately 37 x 22 x 22 Å3, 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 25 independent runs, for which a maximum number of 27,000 generations were produced, employing the Lamarckian Genetic Algorithm (LGA) on the initial population of 50 individuals, a maximum number of 500,000 energy evaluations, with an elitism value of 1, a mutation rate of 0.02, and a cross-over rate of 0.8. 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 1.0 Å 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 4.4.2 [33] 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 (NNB) with a ligand if there is at least one heavy atom of the residue within 3.9 Å of some atom from the ligand. For hydrogen bonds (HHB) the criteria is more restricted. The donor and acceptor atoms of FFRs’ residues and ligands have to be at a maximum distance of 3.5 Å [33].

Automating the molecular docking simulations

We carried out the molecular docking experiments in a Core 2 Quad 2.4 GHz machine, with 8 GB of RAM and 500 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 [33], we also created processing scripts, which are detailed below.

1. For each docked snapshot in a given FFR-ligand complex we extracted and stored the best runs according to the lowest FEB (blind docking) and the lowest RMSD in separate tables. We note that this RMSD measure is calculated with respect to the reference ligand pose (Best Pose in the supervised docking). This produced nine tables (ETH, TCL, and PIF against each of the InhA_wt, InhA_ I16T and InhA_I21V FFR models).

2. 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 36 tables (the previous nine combinations split into HHB and NNB intermolecular contacts).

3. Finally, we counted the amino acid contacts, producing extra 36 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 (SD) of the best FEB (in kcal/mol) and its corresponding RMSD (in Å) with respect to the reference pose for the ligand (Table 2, columns A and B). We also obtained similar statistics for the set of FEB values matching the lowest RMSD with respect to the reference pose for the ligands (Table 2, columns C and D). The latter values are expected to represent the docking results for the best ligand pose (supervised docking). 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 -9.0 ± 2.0, -10.2 ± 1.5, and -10.9 ± 1.4 kcal/mol while for the reference pose, the average FEB was -6.5 ± 2.7, -8.3 ± 2.6 and -8.7 ± 2.3 kcal/mol for InhA_wt, InhA_I16T and InhA_I21V, respectively. The equivalent FEB values for the 1ENY rigid model was -13.4 and -13.5 kcal/mol, respectively. These values represent a difference that varies from 2.5 to 7.0 kcal/mol between the rigid and FFR models of InhA. Also, they are greater than the intrinsic error (2.2 kcal/mol) attributed to the estimation of the FEB by AutoDock3.0.5 [30]. The explicit flexibility of the InhA receptor clearly had an impact in the way PIF interacted with it.

Table 2. Summary of docking results for the FFR models of WT, I16V, and I21T InhA enzyme from M. tuberculosis

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 (Table 3) 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 3 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 80 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 [4].

Table 3. Comparison of the number of different amino acid residues making HHB and NNB contacts with the ligands in at least one of the FFR models’ snapshots during the docking simulations

From the data in Table 3 we were able to single out the top 19 FFR amino acid residues that made contact with ETH and PIF ligands, and the top 18 FFR amino acid residues that made contact with TCL, in at least one of their snapshots (Table 4).

Table 4. Top 19 and 18 amino acid residues that interacted with the ligands in at least one snapshot for each of the three FFR models during the docking simulations

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 (Table 4a) we found three polar (SER20, SER94, THR196), one acidic and one basic (ASP148, LYS165), and 14 hydrophobic residues. As for TCL (Table 4b) we found three polar (THR16, SER94, THR196), one basic (LYS165), and 14 hydrophobic residues. For PIF (Table 4c) we found two polar (SER94, THR196), one acidic and one basic (ASP42, LYS165), and 15 apolar residues.

As can be seen in column 1 of Table 4, most of the residues are hydrophobic, which is somewhat expected since the InhA active site is predominantly apolar. Figure 3 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.

thumbnailFigure 3. Top amino acid residues in the receptor-ligand interactions. Representation of the top 19 and 18 amino acid residues in the receptor structure (PDB ID: 1ENY). The stick representation of ETH (a), TCL (b), and PIF (c) are colored yellow, dark blue, and green, respectively. The residues are represented as magenta beads. The InhA receptor main-chain is represented by ribbons.

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 [7,14,17,34]. We know that the enzyme InhA from M. tuberculosis is a highly flexible receptor [21] and is a very important molecular target for the development of new drugs against tuberculosis [19]. Given that only a few studies [7,12,14,17,34] 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 [21]. 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 (ETH, TCL and PIF) 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 [35], the snapshots at 1.0 ns and 1.5 ns have both 1,559 Å3 and 1,955 Å3, respectively. There is a variation of approximately 400 Å3 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 4 of the nicotinamide portion of NADH to form the ETH-NAD adduct [26], whereas the phenolic ring of TCL forms π-stacking interaction with the NADH nicotinamide ring [27,28]. 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 [26]. Our analyses showed that up to 80 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 80 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 3.1 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 (Bash, Awk and Python). ONS wrote the final version of the manuscript. All authors read and approved the final version of the manuscript.

Acknowledgements

We thank the reviewers for their comments and suggestions which helped to improve the quality of the manuscript and the LAD-PUCRS for CPU time. This project was supported all or in part by grants from CNPq (Processes numbers 554782/2008-1, 302641/2009-2 , 551209/2010-0, 559917/2010-4) to ONS and INCT-TB/CNPq to Prof. Diógenes Santiago Santos and PRONEX 2009/FAPERGS to Prof. Luiz Augusto Basso. EMLC was the recipient of a CNPq MSc scholarship. KSM was supported by a CAPES PhD scholarship. ONS is CNPq Research Fellow.

This article has been published as part of BMC Genomics Volume 12 Supplement 4, 2011: Proceedings of the 6th International Conference of the Brazilian Association for Bioinformatics and Computational Biology (X-meeting 2010). The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2164/12?issue=S4

References

  1. Kuntz ID, Blaney JM, Oatley SJ, Langridge R, Ferrin TE: A geometric approach to macromolecule-ligand interactions.

    J. Mol. Biol 1982, 161(2):269-288. PubMed Abstract | Publisher Full Text OpenURL

  2. Lengauer T, Rarey M: Computational methods for biomolecular docking.

    Curr. Opin. Struct. Biol 1996, 6(3):402-406. PubMed Abstract | Publisher Full Text OpenURL

  3. Koshland DE Jr: Application of a theory of enzyme specificity to protein synthesis.

    Proc. Natl. Acad. Sci.USA 1958, 44(2):98-104. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  4. Koshland DE Jr: The key–lock theory and the induced fit theory.

    Angew. Chem. Int. Ed. Eng 1995, 33(23 24):2375-2378. OpenURL

  5. Sotriffer CA, Flader W, Winger RH, Rode BM, Liedl KR, Varga JM: Automated docking of ligands to antibodies: methods and applications.

    Methods 2000, 20(3):280-291. PubMed Abstract | Publisher Full Text OpenURL

  6. Cozzini P, Kellogg GE, Spyrakis F, Abraham DJ, Costantino G, Emerson A, Fanelli F, Gohlke H, Kuhn LA, Morris GM: Target Flexibility: An Emerging Consideration in Drug Discovery and Design†.

    J. Med. Chem 2008, 51(20):6237-6255. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  7. Carlson HA, McCammon JA: Accommodating protein flexibility in computational drug design.

    Mol. Pharmacol 2000, 57(2):213. PubMed Abstract | Publisher Full Text OpenURL

  8. Erickson JA, Jalaie M, Robertson DH, Lewis RA, Vieth M: Lessons in molecular recognition: the effects of ligand and protein flexibility on molecular docking accuracy.

    J. Med. Chem 2004, 47(1):45-55. PubMed Abstract | Publisher Full Text OpenURL

  9. Cavasotto CN, Abagyan RA: Protein flexibility in ligand docking and virtual screening to protein kinases.

    J. Mol. Biol 2004, 337:209-225. PubMed Abstract | Publisher Full Text OpenURL

  10. Totrov M, Abagyan R: Flexible ligand docking to multiple receptor conformations: a practical alternative.

    Curr. Opin. Struct. Biol 2008, 18(2):178-184. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  11. van Gunsteren WF, Berendsen HJC: Computer simulation of molecular dynamics: Methodology, applications, and perspectives in chemistry.

    Angew. Chem. Int. Ed. Eng 1990, 29(9):992-1023. Publisher Full Text OpenURL

  12. Lin JH, Perryman AL, Schames JR, McCammon JA: Computational drug design accommodating receptor flexibility: the relaxed complex scheme.

    J. Am. Chem. Soc 2002, 124(20):5632-5633. PubMed Abstract | Publisher Full Text OpenURL

  13. Alonso H, Bliznyuk AA, Gready JE: Combining docking and molecular dynamic simulations in drug design.

    Med. Res. Rev 2006, 26(5):531-568. PubMed Abstract | Publisher Full Text OpenURL

  14. Amaro RE, Baron R, McCammon JA: An improved relaxed complex scheme for receptor flexibility in computer-aided drug design.

    J Comput Aided Mol Des 2008, 22(9):693-705. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  15. B-Rao C, Subramanian J, Sharma SD: Managing protein flexibility in docking and its applications.

    Drug Discov Today 2009, 14(7-8):394-400. PubMed Abstract | Publisher Full Text OpenURL

  16. Huang Z, Wong CF: Conformational selection of protein kinase A revealed by flexible ligand flexible protein docking.

    J. Comput. Chem 2009, 30(4):631-644. PubMed Abstract | Publisher Full Text OpenURL

  17. Wong CF: Flexible ligand-flexible protein docking in protein kinase systems.

    Biochim Biophys Acta 2008, 1784(1):244-251. PubMed Abstract | Publisher Full Text OpenURL

  18. Wei BQ, Weaver LH, Ferrari AM, Matthews BW, Shoichet BK: Testing a flexible-receptor docking algorithm in a model binding site.

    J. Mol. Biol 2004, 337(5):1161-1182. PubMed Abstract | Publisher Full Text OpenURL

  19. Agüero F, Al-Lazikani B, Aslett M, Berriman M, Buckner FS, Campbell RK, Carmona S, Carruthers IM, Chan AW, Chen F, Crowther GJ, Doyle MA, Hertz-Fowler C, Hopkins AL, McAllister G, Nwaka S, Overington JP, Pain A, Paolini GV, Pieper U, Ralph SA, Riechers A, Roos DS, Sali A, Shanmugam D, Suzuki T, Van Voorhis WC, Verlinde CLMJ: Genomic-scale prioritization of drug targets: the TDR Targets database.

    Nat. Rev. Drug Discov 2008, 7(11):900-907. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  20. Dessen A, Quemard A, Blanchard JS, Jacobs WR Jr, Sacchettini JC: Crystal structure and function of the isoniazid target of Mycobacterium tuberculosis.

    Science 1995, 267(5204):1638. PubMed Abstract | Publisher Full Text OpenURL

  21. Schroeder EK, Basso LA, Santos DS, Norberto de Souza O: Molecular Dynamics Simulation Studies of the Wild-Type, I21V and I16T of Isoniazid_Resistant Mycobacterium tuberculosis Enoyl Reductase (InhA) in Complex with NADH: Towards the Understanding of NADH-InhA Different Affinities.

    Biophys J 2005, 89:876-884. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  22. Basso LA, Zheng R, Musser JM, Jacobs WR Jr., Blanchard JS: Mechanisms of isoniazid resistance in Mycobacterium tuberculosis: enzymatic characterization of enoyl reductase mutants identified in isoniazid-resistant clinical isolates.

    J. Infect.Dis 1998, 178:769-775. PubMed Abstract | Publisher Full Text OpenURL

  23. Case DA, Cheatham TE III, Darden T, Gohlke H, Luo R, Merz KM Jr, Onufriev A, Simmerling C, Wang B, Woods RJ: The Amber biomolecular simulation programs.

    J. Comput. Chem 2005, 26(16):1668-1688. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  24. Goodsell DS, Morris GM, Olson AJ: Automated docking of flexible ligands: applications of AutoDock.

    J. Mol. Recognit 1996, 9(1):1-5. PubMed Abstract | Publisher Full Text OpenURL

  25. Machado KS, Schroeder EK, Ruiz DD, Norberto de Souza O: Automating molecular docking with explicit receptor flexibility using scientific workflows.

    Lect. Notes Comput. Sc 4643:1-11. OpenURL

  26. Wang F, Langley R, Gulten G, Dover LG, Besra GS, Jacobs WR, Sacchettini JC: Mechanism of thioamide drug action against tuberculosis and leprosy.

    J. Exp. Med 2007, 204(1):73. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  27. McMurry LM, Oethinger M, Levy SB: Triclosan targets lipid synthesis.

    Nature 1998, 394(6693):531-532. PubMed Abstract | Publisher Full Text OpenURL

  28. Kuo MR, Morbidoni HR, Alland D, Sneddon SF, Gourlie BB, Staveski MM, Leonard M, Gregory JS, Janjigian AD, Yee C: Targeting tuberculosis and malaria through inhibition of enoyl reductase.

    J. Biol. Chem 2003, 278(23):20851. PubMed Abstract | Publisher Full Text OpenURL

  29. Oliveira JS, Sousa EH, Basso LA, Palaci M, Dietze R, Santos DS, Moreira IS: An inorganic iron complex that inhibits wild-type and an isoniazid-resistant mutant 2-trans-enoyl-ACP (CoA) reductase from Mycobacterium tuberculosis.

    Chem Commun (Camb) 2004, (3):312-313. PubMed Abstract | Publisher Full Text OpenURL

  30. Oliveira JS, de Sousa EH, de Souza ON, Moreira IS, Santos DS, Basso LA: Slow-onset inhibition of 2-trans-enoyl-ACP (CoA) reductase from Mycobacterium tuberculosis by an inorganic complex.

    Curr Pharm Des 2006, 12(19):2409-2424. PubMed Abstract | Publisher Full Text OpenURL

  31. Humphrey W, Dalke A, Schulten K: VMD: visual molecular dynamics.

    J Mol Graph 1996, 14(1):33-38, 27-28. PubMed Abstract | Publisher Full Text OpenURL

  32. Morris GM, Goodsell DS, Halliday RS, Huey R, Hart WE, Belew RK, Olson AJ: Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function.

    J. Comput. Chem 1998, 19(14):1639-1662. Publisher Full Text OpenURL

  33. Wallace AC, Laskowski RA, Thornton JM: LIGPLOT: a program to generate schematic diagrams of protein-ligand interactions.

    Protein Eng 1995, 8(2):127-134. PubMed Abstract | Publisher Full Text OpenURL

  34. McCammon JA: Target flexibility in molecular recognition.

    Biochim Biophys Acta 2005, 1754(1-2):221-224. PubMed Abstract | Publisher Full Text OpenURL

  35. Dundas J, Ouyang Z, Tseng J, Binkowski A, Turpaz Y, Liang J: CASTp: computed atlas of surface topography of proteins with structural and topographical mapping of functionally annotated residues.

    Nucleic Acids Res 2006, 34:W116-W118. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL