Pharmacophore and Molecular Docking Guided 3D-QSAR Study of Bacterial Enoyl-ACP Reductase (FabI) Inhibitors

Enoyl acyl carrier protein (ACP) reductase (FabI) is a potential target for the development of antibacterial agents. Three-dimensional quantitative structure-activity relationships (3D-QSAR) for substituted formamides series of FabI inhibitors were investigated using comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) techniques. Pharmacophore and molecular docking methods were used for construction of the molecular alignments. A training set of 36 compounds was performed to create the 3D-QSAR models and their external predictivity was proven using a test set of 11 compounds. Graphical interpretation of the results revealed important structural features of the formamides related to the active site of FabI. The results may be exploited for further optimization of the design of new potent FabI inhibitors.


Introduction
The emergence of bacterial resistance to most of the antibiotics currently in clinical use is of world-wide concern [1,2]. In particular, methicillin-resistant Staphylococcus aureus (MRSA) [3] and penicillin-resistant Streptococcus pneumoniae (PRSP) [4] have become troublesome due to the ineffectiveness of remaining therapeutics. Recently, very few novel classes of antibacterial agents have been marketed. Thus, it is urgently desired to develop antibacterial agents with novel mechanism of action against resistant strains.
As fatty acid biosynthesis in pathogenic microorganisms is essential for cell viability, the enzymes involved in the FAS pathway have recently attracted considerable interest as a genomics-driven target for antibacterial drug discovery [5][6][7]. The NADH-dependent enoyl acyl carrier protein reductase (FabI) is a key enzyme in the last step of each cycle of fatty acids elongation [8]. It catalyzes the NADH-dependent stereospecific reduction of α,β-unsaturated fatty acids bound to the acyl carrier protein [9,10]. FabI has been identified to be essential for bacterial viability [8]. In recent years, a wide range of structural classes has been identified as FabI inhibitors [11,12], such as triclosan [13][14][15][16], diazaborines [17,18], imidazoles [19], indole naphthyridinones [20][21][22], thiopyridine [23] and 4-pyridone [24], etc., which demonstrates that FabI is a valid target for antibacterial therapy. Based on indole naphthyridinones, further chemical optimization studies and structure activity relationship studies led to the identification of spiro-naphthyridinone piperidines [25] and pyridodiazepines [26] with improved efficacy and desired physiochemical properties. These inhibitors have shown promising results in the preclinical testing in various resistant strains and mouse models. Despite these and other known inhibitors, more structurally diverse inhibitors of FabI need to be discovered for improving the understanding of the biological function of FabI.
With a view to identify potential compounds with higher predicted potencies, novel FabI inhibitors could be designed by virtual screening and molecular docking based on the available X-ray crystal structures of FabI [19,21,22]. However, unlike quantitative structure-activity relationship (QSAR) techniques, these methods do not readily yield information regarding the importance of the molecular substructures for activity. Comparative molecular field analysis (CoMFA) [27] and comparative molecular similarity indices analysis (CoMSIA) [28] are the most popular three-dimensional quantitative structure-activity relationship (3D-QSAR) methods. Their graphic results can provide a direct way to visualize the structure-activity relationships. To the best of our knowledge, no QSAR studies have been reported by now for those FabI inhibitors. In order to gain further insights into the structural and chemical features required for FabI inhibitory activity, we performed 3D-QSAR analyses of a series of formamides published in the literature [21,22,25,26]. An important challenge for the investigated compounds is to deal with their flexibility, since the applications of CoMFA and CoMSIA require the optimized 3D conformations of all molecules. Alignments of ligands that are known to assume the bioactive conformation were generated by superimposition on pharmacophore points and docking into the protein binding site. The 3D-QSAR results provided useful information about the structural requirements for FabI inhibitors and are expected to usefully support drug design of new FabI inhibitors.

Data Sets
A dataset comprising 47 substituted formamides was taken from the literature [21,22,25,26]. The biological data were considered comparable and divided into a training set and a test set as shown in Table 1. Thirty-six compounds were selected randomly as the training set for model construction and the remaining 11 (asterisked in Table 1) were used as test set for model validation, according to biological activity range and structural diversity. All further studies were performed using the same training and test sets. The IC 50 values were converted into pIC 50 (−logIC 50 ) for use in 3D-QSAR analysis.

Pharmacophore Model Generation
The crystal structure of E. coli FabI with compound 20 (PDB code: 1MFP) was used as starting structure for the generation of the pharmacophore model. The software LigandScout 3.01 [29,30] was used for detection and interpretation of crucial interaction patterns between FabI and the ligand. LigandScout extracts and interprets ligands and their macromolecular environment from PDB files and automatically creates and visualizes an advanced pharmacophore model. Then the pharmacophore model was exported as a hypoedit script and converted into Discovery Studio 2.1 [31] format with Hypoedit tool. Subsequently, the pharmacophore model was used for mapping all of the molecules.

Molecular Docking
The docking procedure aims to generate and score putative protein-ligand complexes according to their calculated binding affinities. Docking studies were carried out using GOLD docking software [32], version 3.1, which uses a powerful genetic algorithm (GA) method for conformation search and docking, and is widely regarded as one of the best docking programs [33]. Docking experiments were performed using the default GOLD fitness function (VDW = 4.0, H-bonding = 2.5) and default evolutionary parameters: population size = 100; selection pressure = 1.1; operations = 100,000; islands = 5; niche size = 2; migration = 10; mutation = 95; crossover = 95. The ChemScore function was used to rank different binding poses. The center of the bound ligand was defined as the binding site. Ten docking runs were performed per structure. All poses were output into a single *.sdf file.

Alignment Rule
In the 3D-QSAR studies, the molecular alignment and conformation determination are very important to construct reliable models. Due to the flexibility of the investigated compounds, it is difficult to choose a suitable conformation that achieves a meaningful superimposition. In an ideal alignment the biologically active conformations should be aligned taking into account the orientations that the ligands adopt at the binding site of the protein. Therefore, we applied two different receptor-based alignments, with the conformations obtained from structure-based pharmacophore (SBP) search and docking.
All the molecules in the training and test sets were mapped simultaneously onto the pharmacophore model using the "flexible" fitting method and "best mapping only" option in the Ligand Pharmacophore Mapping protocol in Discovery Studio 2.1. The conformation selected for each compound, which was assumed to be the bioactive conformation, corresponded to the conformation which best fit the pharmacophore model. The final aligned molecules were exported to SYBYL6.9 [34] for CoMFA and CoMSIA analysis.
For the docking, all the molecules were docked into the FabI active site using the GOLD program. The conformation with the highest ChemScore of each molecule and their alignment were used directly in CoMFA and CoMSIA to explore 3D-QSAR models.

CoMFA and CoMSIA Model
CoMFA was performed using the QSAR option of SYBYL 6.9. The steric and electrostatic field energies were calculated using the Lennard-Jones and the Coulomb potentials, respectively, with a 1/r distance-dependent dielectric constant in all intersections of a regularly spaced (0.2 nm) grid. The electrostatic fields were computed using Gasteiger-Huckel charge calculation methods. A sp3 hybridized carbon atom with a radius of 1.53 Å and a charge of +1.0 was used as a probe to calculate the steric and electrostatic energies between the probe and the molecules using the Tripos force field. The standard parameters implemented in SYBYL 6.9 were used. The truncation for both the steric and the electrostatic energies were set to 30 kcal/mol.
The CoMSIA method involves a common probe atom and similarity indices calculated at regularly spaced grid intervals for the pre-aligned molecules. The similarity indices descriptors were derived with the same lattice box used in CoMFA. CoMSIA calculates hydrophobic, H-bond donor and acceptor fields in addition to steric and electrostatic fields. The distance dependence between the grid point and each atom of molecule was determined by the Gaussian function through the similarity indices calculated at all grid points, and a default value of 0.3 was used as an attenuation factor.

Partial Least Square Analysis
Partial least squares (PLS) [35] methodology was used for all the 3D-QSAR analyses. The cross-validation [36,37] analysis was performed using the leave one out (LOO) method in which one compound is removed from the dataset and its activity is predicted using the model derived from the rest of the dataset. The cross validated r 2 that resulted in the optimum number of components and lowest standard error of prediction were considered for further analysis. To speed up the analysis and reduce noise, a minimum filter value σ of 2.00 kcal/mol was used. Final analysis was performed to calculate conventional r 2 using the optimum number of components obtained from the cross-validation analysis.
The predictive power of the 3D-QSAR models was determined from a set of eight molecules that were excluded during model development. The optimization, alignment and all other steps of these test set molecules were the same as those of the training set molecules described above, and their activities were predicted using the model produced by the training set. The predictive correlation (r pred 2 ) based on the test set molecules, is computed using where, SD is defined as the sum of the squared deviations between the biological activities of the test set and mean activity of the training set molecules and PRESS is the sum of the squared deviation between the predicted and actual activity values for each molecule in the test set.

Pharmacophore Elucidation
As shown in Figure 1a, the pharmacophore model automatically generated by the LigandScout 3.01 program includes five features: three hydrogen bond acceptors (HBA), one hydrogen bond donor (HBD) and one hydrophobic group. Besides, the program automatically generated a series of excluded volumes in the model. Two HBA features characterize the carbonyl group of the ligand, which forms two hydrogen bonds with Tyr156 and the 2'-hydroxyl group of the nicotinamide ribose of the nucleotide (Figure 1b). The other HBA and the HBD features are involved in H-bond interactions between alanine 95 and the pyridylamine and N-acyl hydrogen of the naphthyridinone functionality, respectively. The hydrophobic feature is located on the indole portion. To verify whether the pharmacophore model finds the correct bioactive conformation, we applied the method to two structurally similar compounds 20 and 9, whose bioactive conformations are known from E. coli X-ray structures of the ligand-enzyme complex. Their bound conformations were mapped onto the pharmacophore model using the "flexible" fitting method and "best mapping only" option in the Ligand Pharmacophore Mapping protocol and then superimposed on the best mapping conformations (Figure 2). From the results, it is important to note that the best mapping conformations of the ligands are in agreement with their bioactive conformations. Hence, the results showed that the pharmacophore model is capable of reproducing the bioactive conformation from the Protein Data Bank and is reliable enough to retrieve compounds that fit all the features of the query from chemical databases. The final aligned molecules are shown in Figure 3a.

CoMFA and CoMSIA Results
The stepwise development of CoMFA and CoMSIA models using different fields is presented in  The external validation is important to establish a reliable 3D-QSAR model. To validate the stability and predictivity of the 3D-QSAR models, 11 compounds that were not included in the construction of the 3D-QSAR models were selected as the test set. The CoMSIA-SEHD model showed the best reasonable external predictivity, yielding a r pred 2 of 0.702 ( Table 2). The predicted values of training and test sets are presented in Table 3. The correlations between the predicted and experimental values of all compounds are shown in Figure 4a.

Docking Analysis
To test whether the GOLD program is feasible for ligand binding to FabI, two X-ray crystal structures (PDB code, 1FMP [22] and 1LXC [21]) were initially chosen and the docked conformation of ligands was compared with their crystallographic conformation. The proteins were prepared by removing the small ligands and adding hydrogens using SYBYL 6.9. The active site for docking was defined as all atoms within 8 Å radius of the co-crystallized ligand. The resulting docked conformation with highest ChemScore value and that of crystallographic conformation were very similar ( Figure 2). This result demonstrates that GOLD analysis is suitable for the identification of the binding mode between inhibitors and FabI. The resulting alignment of all compounds is shown in Figure 3b. This molecular alignment is more realistic and more suitable for molecular field analyses.

CoMFA and CoMSIA Results
The CoMFA and CoMSIA results generated from docked conformations are presented in Table 4. The CoMFA models showed r cv 2 of 0.664 with six components, r 2 of 0.993 and F value of 730.83. For CoMSIA, differing from the pharmacophore alignment, the best r cv 2 was obtained using only steric and electrostatic descriptor variables. The high value of r cv 2 appears to be a necessary but not the sufficient condition for the model to have a high predictive power. A CoMSIA-SEHD model with a r cv 2 of 0.711 with six components, a r 2 of 0.995 and a best external predictivity r pred 2 of 0.864 was obtained. These data demonstrate that a more statistically robust model was obtained from the CoMSIA study. The experimental versus predicted activities (Table 4) using CoMSIA-SEHD models of all compounds are shown in Figure 4b.

CoMSIA Models Interpretation
Steric, electrostatic, hydrophobic and hydrogen bond donor contour plots obtained using CoMSIA analysis from pharmacophore and docking based alignment are useful to explore the protein-ligand interactions, presented in Figures 5 and 6, respectively. For simplicity, the interactions between only compound 20 and the contour plots are shown. As shown in Figure 5a, the green contour observed above the naphthyridinone ring indicated that some bulky substitutions at this region would be favorable for the activity, which was known to be exposed to solvent. This can explain the fact that compounds 26-33 with spiropiperidine substitutions and compounds 34-47 with seven-and eight-membered heterocyclic rings exhibited potent FabI inhibitory activity. A large yellow contour found around the ene-amide of compound 20 suggested that any bulky substitution at this position would be likely to decrease activity. The orientation of the linking amide may be very important for activity.
The electrostatic contour map is shown as red and blue contours in Figure 5b. The small blue contour near the naphthyridinone ring suggested that the presence of positively charged groups in this region was favorable for activity. Besides, a large positively charged blue contour was found near ene-amide too. Two carbonyls of compound 20 were oriented towards two red contours, thereby indicating favorable interactions with negatively charged groups.
The hydrophobic contour map of CoMSIA is presented in Figure 5c. The hydrophilic favored white contour around the naphthyridinone ring implied that any hydrophilic substitution was favored there.
This was in accordance with compounds 25-33, where incorporating hydrophilic morpholine or piperidine proved to improve the activity. The yellow contour indicated that any hydrophobic group substituent here was favored. This region supported the observation that compounds 29 and 33, with OPr substitution, were among the compounds with a high FabI inhibitory activity.
The hydrogen bond donor contour plot is displayed in Figure 5d. The cyan contours represent the region favoring hydrogen bond donor substituents. The purple contours indicate that hydrogen bond acceptor groups in these regions are required for higher activity. It can be explained by the fact that the hydrogen bonds were involved in the naphthyridinone functionality and enzyme FabI. That was why compounds 1-16 without naphthyridinone were less potent compared to compounds 19-24, etc. For compounds 17 and 18, acylated amino pyridine was less favored compared to the hydrogen bond of naphthyridinone.

Molecular Docking-Based Alignment
The CoMSIA contour plots from the docking-based alignment are shown in Figure 6. The steric, electrostatic and hydrophobic contour maps of CoMSIA were similar to those of CoMSIA from the pharmacophore-based alignment in explaining the effect of substitution on the biological activity (Figure 6a-c). One difference was that a yellow region around the pyridine ring represented an area where no substitution was favored.
For the hydrogen bond donor contour map, the two small cycan contours around the N-acylhydrogen indicated that hydrogen bond donor substituents would improve the activity. This can be explained by the fact that compounds 17-47, which possessed N-acylhydrogen, showed better potencies than compounds 1-16. The other cycan contour above the naphthyridinone revealed that hydrogen bond donor substituents might be favored. This may be the reason why a series of seven-and eight-membered heterocyclic derivatives 34-39 and 41-47 bearing a free NH atom showed improved potencies compared to compound 40 with nitrogen methylation. A small purple contour at the nitrogen position of pyridine suggested that a hydrogen bond acceptor substituent at this position might be favored. Most of the active derivatives (11, 12, 14-47) possessed hydrogen bond acceptor groups such as a nitrogen atom at this position. Compounds 7-10 and 13 without a hydrogen bond acceptor substituent at this position exhibited significantly decreased potencies. Another purple contour around the naphthyridinone also revealed that a hydrogen bond acceptor group would improve the activity. This was consistent with the fact that compounds 45 and 46 bearing the carbonyl group at this position exhibited a higher activity than compound 47. Despite the different alignments, their 3D-QSAR models are comparable to those presented there. However, this docking-based study provides significantly better predicted statistical coefficients (r pred 2 = 0.864) and additional information about the hydrogen bond interactions of the naphthyridinone moiety, which emphasizes the benefit of combining docking for the selection of the alignment followed by 3D-QSAR for formamides-based FabI inhibitors.

Summary of Structure-Activity Relationships
Comparing the 3D-QSAR models, the improved performance of the CoMSIA is obvious, while the differences in the case of CoMFA are negligible. The results of 3D-QSAR CoMSIA obtained from both alignments are complimentary to each other. The detailed contour analysis of both CoMSIA models enabled us to point out several structural requirements as mentioned in the above discussion for the observed inhibitory activities (Figure 7). In detail, the bulky and hydrophilic groups at the S 1 region are favorable, while the hydrophobic groups at the S 2 region are favorable. The hydrogen bond acceptor group at the H 1 and H 3 position may benefit the potency. The hydrogen bond donor group at the H 2 position increases the activity. The ene-formamide skeleton was essential for high activity.

Conclusions
In this study, the conformations of all compounds obtained from the pharmacophore-based alignment and docking into the active site of FabI were used for the CoMFA and CoMSIA analyses. Both alignment procedures led to statistically robust and predictive 3D-QSAR models. Despite the different alignments, the contour maps from both models are similar in explaining the influence of substitution on activity. The predictive ability of the models was confirmed by predicting the activity of 11 compounds used as external test sets. The robust 3D-QSAR model and its three-dimensional contour map provide guidelines to design compounds with new scaffolds and optimize known molecules.