In Silico Studies of Quinoxaline-2-Carboxamide 1,4-di-N-Oxide Derivatives as Antimycobacterial Agents

Molecular modelling studies were performed on some previously reported novel quinoxaline-2-carboxamide 1,4-di-N-oxide derivatives (series 1–9). Using the LigandScout program, a pharmacophore model was developed to further optimize the antimycobacterial activity of this series of compounds. Using the Dock6 program, docking studies were performed in order to investigate the mode of binding of these compounds. The molecular modeling study allowed us to confirm the preferential binding mode of these quinoxaline-2-carboxamide 1,4-di-N-oxide derivatives inside the active site. The obtained binding mode was as same as that of the novobiocin X-ray structure.


Introduction
Tuberculosis (TB) is the most prevalent infectious bacterial disease caused by Mycobacterium tuberculosis (mtb). In 2009, the World Health Organization (WHO), reported an estimated 9.27 million OPEN ACCESS cases of mtb in 2007. A remarked increase was noticed from the 9.24 million cases in 2006, the 8.3 million cases in 2000 and the 6.6 million cases in 1990. While the total number of incident cases of mtb ius increasing, the number of mtb infected cases per capita is slowly decreasing. A leading killer, tuberculosis is an intracellular infection responsible for some 3 million deaths annually, with a person lost to mtb every 15 s [1]. Readily spreading from person to person, and showing bad resistance to isoniazide and rifamicin, multidrug-resistant strains of M. tuberculosis (MDRTB) will necessarily make the future control of TB more difficult. This dilemma was worsened by the emergence of XDR-TB, that is unresponsive not only to isoniazid and rifampicin, the first line drugs, but also to a fluoroquinolone and to at least one of the second-line drugs (amikacin, capreomycin or kanamycin) [2]. In addition to these drug-resistance characteristics of TB, the recent influx of immigrants from countries endemic for the disease and co-infection with human immunodeficiency virus (HIV) [3,4] highlight the urgent need for new drugs to extend the range of effective TB treatment options. Quinoxaline derivatives are an attractive class of target compounds for new drug development because of their potentially versatile biological activities, including antiviral, anticancer, antibacterial and antiprotozoal properties [5][6][7][8][9]. As antituberculous agents, a wide range of quinoxaline-1,4-di-N-oxide derivatives with variable substituents at different positions were reported [10][11][12][13][14][15][16][17][18].
As anti-T. cruzi agents, QSAR studies of 3-arylquinoxaline-2-carbonitrile di-N-oxides were performed by establishing a link between IC 50 values and their Moriguchi-octanol/water partition coefficients (MLOGPs) [19]. In the past decades, due to a lack of experimental biological data on quinoxaline derivatives as anti-mycobacterial agents, none of the studies have included quantitative structure activity results. With the addition of new data to the literature, research groups have undertaken further exploration of the biological profile exhibited by quinoxaline-2-carboxylate 1,4-di-N-oxide derivatives resorting to the QSAR formalism, establishing predictive models for biological properties [20].
In view of these findings, we are prompted to perform further studies, using molecular modeling software, in order to explore the structural requirements necessary for the anti-tuberculous activity of quinoxaline derivatives. Furthermore, in silico interactions of the quinoxaline derivatives under study within the mtb-DNA gyrase active site could provide valuable information for their possible mode of action.

Pharmacophore Modeling
Elucidation of the binding approaches for the compounds under study is suggested based on finding the active structures. Figure 1 shows the structure of the training set compounds (1b,c,f, 2a,b,d,g, 3a, 4b-e, 5a,e, 6a,b, 7b,e,g, 9a,b,e) and Figure 2 shows the structure of the test set compounds (1a,d,e,g, 2c,e,f, 3b, 4a,f,g, 5b,g, 6e,g, 7a, 8a,b,e,g, 9g). Based on the assumption that the active compounds bind in a similar fashion at the active site, the Ligandscout program [19][20][21] was employed to evaluate the common features essential for antiproliferative activity and the hypothetical geometries adopted by these ligands in their most active forms. Thus, these compounds were submitted for pharmacophore model generation based on the shared chemical features. Diverse conformations within a 20 kcal/mol energy range were generated and submitted to the alignment procedure.  Structures of the quinoxaline-2-carboxamide 1,4-di-N-oxide compounds used in the LigandScout test set [16].
The successful pharmacophore run resulted in generation of 10 hypotheses, as its highest rank score and mapping into all training set molecules, hypo1 was considered statistically as the best hypothesis and it was selected for further investigation and analysis. The top-ranked chemical feature-based pharmacophore model identified in this study is shown in Figure 3. This pharmacophore model contains nine chemical features: one aromatic ring (blue), four hydrophobes (yellow), three hydrogen acceptors (red) and one hydrogen donor (green). All the training set and test set compounds were mapped onto hypo1 with scoring the orientation of a mapped compound within the hypothesis features using a "fit value" score. Figure 3. The top-ranked chemical feature-based pharmacophore model developed using the LigandScout program. The pharmacophore includes one aromatic rings (blue), four hydrophobes (yellow), three hydrogen acceptor feature (red) and one hydrogen donor (green). Distances are given in Angstrom.
As a quick and primary validation of hypo1, mapping of the compounds found to show a good agreement between the fit value and the biological activity (Tables 1 and 2 Tables 1 and 2 revealed a moderate correlation between the fit value and the biological activity of each of the tested compounds. Of the training set, the highly active compounds (1b, 2b, 4b-d, 6b and 9b) showed a range of fit value of 131.46-132.89 whereas compounds with lower activity showed a lower fit value average of 125.6-114.85. This initial correlation encouraged us to generate a linear model based on "fit value" to predict the biological activity of the compounds under investigation. The generated model (Equation 1) showed good statistics and was used successfully to calculate the activity of the tested compounds ( Figure 2):

, Figures 4 and 5). Initial investigation of the results shown in
where number of compounds, x is the fit value and n = 22, st. error = 0.1895, R= 0.929 and R ² = 0.863.     (Figure 6a). A slight displacement of the phenyl ring substituent on the amide moiety away from the plane of the aromatic pharmacophore center results in variation of activity (Figure 6b). At the C-7 of the quinoxaline ring, the displacement of the hydrophobic substituent away from the hydrophobic pharmacophore center (Figure 6b) or the presence of non-hydrophobic groups (Figure 6c) can partially explain their lack of activity. The rest of the features that are common for all compounds are the oxygen atoms adjacent to N1 and N4 positions of the quinoxaline scaffold, as hydrogen acceptor, the amide NH as hydrogen bond donor and the methyl or phenyl group at C-3 position.

Homology Modeling
Homology modeling and subsequently, a docking process were undertaken in order to inspect the prospective interactions between the quinoxaline 1,4-N-dioxide derivatives, and the active site of the Mycobacterium DNA gyrase B subunit. In the homology modeling study, the crystal structure of the gyrase B 43 K ATPase domain complex with the potent inhibitor novobiocin (1KijB.pdb) [22], was selected as the template structure. This particular template has been selected based not only on BLAST-p alignment but also on the structural similarity between our quinoxaline 1,4-di-N-oxide derivatives and co-crystallized novobiocin. During the homology modeling process the protein coordinates were first minimized using the AMBER94 force field, then the heavy atoms were modeled and followed by addition of all hydrogen atoms. The pair-wise percentage residue identity was determined as 41.321% between two chains, where the pair-wise RMSD values, for Cα atoms of the superimposed model and template, was 0.611 Å. In brief, the model structure comprises a compact single domain including an 8-stranded beta sheet and 6-alpha-helices and random coils (Figure 7).

Docking Procedure
A docking study was undertaken using Dock6.4 [23] in order to investigate the possible interactions between the designed compounds and the active site of the epidermal growth factor receptor (EGFR) and to compare it with the binding mode of the known 1kijB inhibitor, novobiocin. The X-ray structure of the enzyme bounded with novobiocin was taken from the protein data bank (PDB code: 1kijB) [22]. The RMSD value difference of 0.671415 Å of the pose of the non-restricted redocking of the X-ray structure of the gyrase inhibitor (novobiocin) from itself also confirmed the approach (Figure 8). The docking poses of compound 4d, as an example of the designed compounds compared with that of the novobiocin, X-rary structure, are shown in Figures 9-11, respectively. Occupying the same binding site as well as the coumarin moiety of the novobiocin X-ray structure, the quinoxaline-1,4-dioxide scaffold structure showed a different orientation. The substituted benzyl ring of the amide moiety is oriented in a binding pocket surrounded by ASN8, Ala9, Ala18, Gly27, ARG28, Gly29, Ile30, Gly83, Ser84 with hydrogen bonding between the oxygen atom of the amide group and the amino group of the side chain of Arg77.

General
All molecular modeling studies were performed on PC Windows Vista Home Premium Intel(R) Core(TM)2 Duo, 1.83 GHz using the LigandScout program v3.1 1999-2013 (G. Wolber and Inte: Ligand GmbH, Vienna, Austria) [21] and the Dock6.4 program [23]. The quinoxaline 1,4-di-N-oxide derivatives used in this study are shown in Figures 1 and 2. The biological data are cited from the literature [16].

Ligand Based Pharmacophore Modelling
The study was carried out using the software LigandScout (version 3.0). Using the default settings, the LigandScout program [21] was used to derive the 3D chemical feature-based pharmacophores from the structural data of the compounds (1, 2, 4) a-g, 3a,b and (5-9) a, b, e, g [16] included in the modeling method. Prior to the generation of the pharmacophore hypotheses, the training set  compounds 1b,c,f, 2a,b,d,g, 3a, 4b-e, 5a,e, 6a,b, 7b,e,g, 9a,b,e (Figure 1) were converted to 3D structures and used to generate diverse conformations. The diverse conformation generation protocol implemented in the LigandScout program was used to generate conformations using the best conformation model generation method. Other parameters like maximum number of 500 conformers, and an energy threshold value of 20 kcal/mol above the global energy minimum were chosen during conformation generation. During pharmacophore hypothesis generation four pharmacophoric features like hydrogen bond acceptor (HBA), hydrogen bond donor (HBD), ring aromatic (RA) and hydrophobic (HY) were selected based on the feature mapping results. All parameters were set to their default values.

Pharmacophore Validation
The generated pharmacophore hypothesis was validated using leave-one-out and test set methods.

Leave-One-Out Method
The pharmacophore hypothesis is cross validated by leave-one-out method. In this method, one compound is left in the generation of a new pharmacophore model and its affinity is predicted using that new model. The model building and estimation cycle were repeated until each compound was left out once [24]. This test was performed to verify whether the correlation coefficient of the training set compounds is strongly depend on one particular compound or not [24].

Pharmacophore Validation using Test Set
Compounds 1a,d,e,g, 2c,e,f, 3b, 4a,f,g, 5b,g, 6e,g, 7a, 8a,b,e,g, 9g ( Figure 2) were selected as a test set. This method is used to elucidate whether the generated pharmacophore hypothesis is proficient at predicting the activities of compounds other than the training set and classifying them correctly in their activity scale. The conformation generation for the test set compounds was carried out in a similar way, like the training set compounds using BEST conformation analysis algorithm, implemented within the LigandScout program with setting values, as same as those used with the training set. The compounds associated with their conformations were subsequently carried out for pharmacophore mapping [24].

Homology Modeling
The building of the binding cavity in the specific and original size based on 1kijB (16-427) [22] was successfully obtained using BLAST alignment of the sequence of M. tuberculosis DNA gyrase subunit B (access number P0C5C5) [25] with the 1kijB.pdb sequence structure and followed by HH Search.
Believed to be responsible for key interactions, the crystallized novobiocin structure and a water molecule were kept in their original positions for rebuilding during the modeling process. The BLAST module parameters were selected to run HHSEARCH of 50, minimal number of uncovered target residues to model an additional template of 25 and automated mode of SMR-pipeline. The obtained homology modeled-structure of mtb-DNA gyrase subunit B was used for the docking study in preparation of the input receptor files within Dock6.4 program.

Conclusions
In conclusion, a computational approach along with the 3D-QSAR and docking analysis was employed to identify molecular structural features required for effective antimycobacterial activity, with the aim of discovering drugs for treatment of M. tuberculosis infection. A reliable pharmacophore model was generated based on 22 training set compounds, which consists of one aromatic ring (A), four hydrophobes (HB), three hydrogen acceptors (HA) and one hydrogen donor (HD). This model revealed internal (R 2 = 0.863) prediction of training set as well as external (Predr2 = 0.74) prediction of 21 compounds of test set. The quinoxaline-1,4-N-dioxide scaffold of the docked compound 4d occupied the same binding site as the coumarin moiety in the novobiocin X-ray structure. Compound 4d and the X-ray structure showed similar hydrogen bonding to the amino group of the side chain of Arg77. These findings could be exploited for future ligand design in order to obtain novel derivatives as inhibitors of Mycobacterium tuberculosis.