Hybrid Pharmacophore Approach as a Novel Strategy to identify Anti-Mycobacterial Chemotypes: Investigation with DapB Enzyme

The current study examines the efficacy of dynamics-based hybrid pharmacophore models (DHPM) based on newly explored interaction features, as tools to screen potential inhibitors of Mycobacterium tuberculosis (Mtb)-DapB, a validated target essential for L-Lysin biosynthesis. Molecular dynamics (MD) simulations were performed on Mtb-DapB models, generated from a reported crystal structure by linking the cofactor (NADH) and substrate mimetic inhibitor (2,6-PDC), thereby creating a hybrid molecule (HM). Comparative investigation of the dynamics of interactions made by NADH, 2,6-PDC and HM from the MD trajectories revealed a number of stable interactions between HM and the hinge region residues leading to significant alterations of the Mtb-DapB structure. These newly identified interactions, translated into DHPM may be utilized as Mtb-DapB specific screening filters and provide new insights for structure activity relationships. To validate the applicability of the DHPM, about 10 million compounds compiled from the publicly available chemical space were screened using the DHPM as well as conventional pharmacophore models (based on the NADH/2,6-PDC). Systematic analyses of structures, physicochemical/ADMET properties and molecular docking of the 387 and 982 compounds screened by the DHPM and conventional models respectively demonstrate the capabilities of DHPM to screen structurally diverse chemical entities with better druglike properties and higher target binding potentials. The study demonstrates application of DHPM to predict 8 highly active anti-Mtb (MIC<=2{micro}M) compounds whose mechanisms of action were previously unknown.


Abstract
The current study examines the efficacy of dynamics-based hybrid pharmacophore models (DHPM) based on newly explored interaction features, as tools to screen compounds whose mechanisms of action were previously unknown.

Introduction
Tuberculosis (TB) is the leading cause of death worldwide due to a single infectious agent according to the latest world health organization (WHO) reports [1,2]. The standard therapies of treating TB with a combination of several lines of antibiotics over a period of six to nine months [3,4] have become ineffective due to emergence of drug resistant Mycobacterium tuberculosis (Mtb) strains [5]. The "End Tuberculosis Strategy" [1] of WHO calls for intensified research and innovation in TB drug discovery using multidisciplinary approaches to identify novel drug targets as well as fast and accurate techniques to design new chemical entities with higher potency to address the scourge of Mtb. Computational methods have become essential in last few decades not only to understand the drug-target interactions, but also for in silico screening of huge chemical libraries providing a fast and less expensive alternative to the traditional high throughput screening [6][7][8][9][10][11]. Recent literature has reported the identification of novel Mtb drug targets like DprE1, MmpL3, ATP, mycolic acid and Diaminopimelate (DAP) biosynthesis pathway enzymes etc that may be inhibited to tackle MDR and XDR strains [4,12]. The DapB enzyme of the DAP biosynthetic pathway from H37Rv strain of Mtb (Mtb-DapB) is one of the validated drug targets [13,14] as inhibition of this enzyme blocks the production of meso-diaminopimelate thus leading to inhibition of de novo lysine biosynthesis and peptidoglycan assembly, both of which are crucial for the survival of the pathogen [15] (Fig 1a).
Several groups made efforts to identify inhibitors of enzyme Mtb-DapB by exploring the potential of product analogues [16] as potential inhibitors, but have met with very limited success. 2, 6-PDC (Pyridine-2, 6-dicaboxylate) and few other heterocyclic aromatic product analogues have been identified with IC 50 as high as 26 µM for Mtb-DapB [16]. Also, a number of sulphonamide inhibitors of Mtb-DapB which inhibited Mtb-DapB competitively with respect to the substrate 2,3dihydrodipicolinic acid were identified with K i values ranging from 7 to 48 µM [17].
This indicates that only substrate or product mimicry is not sufficient and new lead generation strategies should be employed exploring the key features of the binding sites of the enzymes. Molecular dynamics (MD) based pharmacophore models have emerged as quite powerful tools which not only account for the flexibility of the targets, but also help to identify novel key chemical features in the binding sites which is otherwise unexplored in the crystal structures and may be harnessed to design new inhibitors [18,6,7].
In this study, we have investigated dynamics-based hybrid pharmacophore modeling as a novel strategy to identify new anti-TB chemical space. MD simulations were performed on two model systems of Mtb-DapB to attain a detailed picture of the non-covalent interactions in the binding site responsible for stable cofactor and substrate binding and explore new regions of the binding site that may be targeted for effective enzyme inhibition. MD based e-pharmacophore models were generated representing the stable interactions of 2,6-PDC (a mimetic of the natural substrate of Mtb-DapB) and NADH (natural cofactor of Mtb-DapB) individually with their respective binding sites as well as hybrid e-pharmacophore models which represented interactions of both 2,6 PDC and NADH. A comprehensive dataset of above 10 million compounds compiled from publicly available chemical space was screened using these models. The structural and physicochemical properties of compounds screened by the hybrid models and the substrate/cofactor-based models were compared to examine the abilities of the hybrid models to identify a newer chemical space with better druglike properties and higher binding affinities with the target.
Although in this article, we have applied the dynamics-based hybrid pharmacophore strategy for Mtb-DapB, it can also be tested with other drug targets.

Model Systems and Molecular Dynamics (MD) Simulations
Two different model systems of Mtb-DapB were considered for our study. The initial coordinates for the first model system (DapB-N-P) were taken from the crystal structure 1P9L [19], binding NADH and 2,6-PDC. The second model system (DapB-Hyb) was generated by connecting the nicotinamide part of NADH with 2,6-PDC of 1P9L using a propyl linker (Fig 1b) thereby creating a hybrid molecule (HM). The two systems were subjected to MD simulations using AMBER12 package [20]. The protein parts were prepared implementing AMBER ff99SB [21] force field while the ligands were prepared using General AMBER Force Field (GAFF) [22]. The partial atomic charges for NADH and PDC were derived using the antechamber module [23] by semi-empirical AM1 method with bond charge correction (BCC) model. A 1000step steepest descent (SD) minimization was carried out on both the system to get rid of nonphysical contacts. Both the systems were then solvated in a cubical TIP3P water box and were neutralized by adding appropriate number of counter ions. A gradual heating from 0 to 300K was performed for 50 ps followed by 500 steps of SD and 500 steps of conjugate gradient energy minimizations. The systems were subjected to constant pressure equilibration (300K and 1 atm pressure) of 40 ns.
Finally, the production run of 40 ns was performed under NPT ensemble with leapfrog integrator and a time step of 0.002 ps. Langevin piston algorithm, SHAKE, and particle mesh Ewald (PME) were used to perform pressure control, constraint covalent bonds involving hydrogen atoms and to treat long-range interactions, respectively. The coordinates were saved every 2 ps for analyses.

Generation of e-pharmacophore models
Eight snapshots were extracted at every 5 ns interval from each of the MD trajectories of the two model systems of Mtb-DapB, which were further used to construct the structure based e-pharmacophore models. Glide energy grids were generated for each snapshot to define the active site as a cubical box of 12 Å × 12 Å × 12 Å around the ligands (NADH and 2,6-PDC in case of DapB-N-P and HM in case of DapB-Hyb) and their interactions with the receptor (Mtb-DapB) were evaluated by using "Score in place" mode of the Glide module of S c h r o dinger molecular modeling suite [24,25]. Default settings, with the option to output Glide extra precision (XP) [25] descriptor information, were used for the scoring. The resulting protein−ligand complexes along with the XP energy terms for hydrophobic enclosure, hydrophobically packed correlated hydrogen bonds, electrostatic rewards, π -π stacking, cation-π, and other interactions were then submitted to the Phase [26] module of Schrodinger to generate energy-based pharmacophore (e-Pharmacophore) models [27]. Each interaction is represented by a pharmacophore feature site and is assigned an energetic value equal to the sum of the Glide XP contributions of the atoms comprising the site. Then, the sites are ranked based on the energetic terms.  Fig 1) were used for H type e-pharmacophore model generation. screening study. This huge chemical space was preliminarily screened by applying the first level ligand based ElectroShape filters [28] implemented in the SwissSimilarity search server [29] taking NADH, 2,6-PDC and HM* as query molecules. A unique set of 6781 molecules screened by ElectroShape method using all the three query molecules were further screened by the dynamics structure-based e-Pharmacophore models. All these 6781 compounds were minimized using the default parameters of LigPrep [30] module of the Schrodinger Suite and five lowest energy conformers were retained for each compound. Screening was performed using the "Advanced Pharmacophore Screening" option of Phase module of Schrodinger Suite [26]. Five conformations per rotatable bond were generated for each ligand and maximum number of conformations per compound was assigned to be 100. A thorough sampling was used for screening and the default option for skipping structures with more than 15 rotatable bonds was used. The minimum number of sites the molecule must match was assigned to be 5 for the N-Type and H-type models, while it was 2 the P-type models. Among many conformers of a ligand, the one with the best fitness score (S) evaluated by a specific fitness function (S1 Equation) [27] was retained for each compound. The compounds screened by the N-, P-and H-type models were named as N-set, P-set and H-set respectively. Various cheminformatics analyses were performed to compare the structural, physicochemical and ADMET properties of the three sets of molecules.

Docking
The N-set and H-set molecules were docked to the MD snapshots obtained from DapB-N-P and DapB-Hyb respectively taking the same grids generated for e-

Screening of Mtb active (MA) compound library
A library of 2950 antitubercular compounds was carefully curated from literature [31] and ChEMBL database which is named as Mtb active (MA) compound library. These compounds have reported low micro molar MICs (Minimum inhibitory concentration) as low as 2 µM in the whole cell assays on Mtb. These compounds were prepared and screened against the eight H-type models with using the same parameters as described in the 'Advanced pharmacophore screening' section. The compounds screened by each H-type model (compounds that matched at least 5 pharmacophore features) were docked to the respective DapB-Hyb snapshots (grids generated for e-pharmacophore modelling were used) with the same Glide XP parameters as described in the 'Docking' section.

Overall Structure and dynamics of Mtb-DapB in DapB-N-P and DapB-Hyb
The Mtb-DapB is a homo-tetramer of 245-residue monomers (Fig 3a). Two major domains and two short hinge regions connecting them comprise each monomer.
The N-terminal domain (1-106 and 216-245) contains six β -strands (β1-β5 and β 10) and four α -helices (α1-α3, and L10(212-215) connect two domains and act as hinge regions for the domain movements. The ADP part of NADH is embedded in a solvent exposed groove like region located in the N-terminal domain and extending to the hinge region, while the nicotinamide part is placed in the floor of a relatively less exposed cavity C1 (Fig 3b) formed by residues from both N-and C   (Fig S1b), which suggests that the individual systems do not show major structural deviations during the simulations and are well equilibrated. Stable energy profiles of the systems also support the above observation (Fig S1c). Hence, the last 10 ns of both the trajectories were used for further analyses. The plots showing the interaction energies between the ligands and receptor (Fig S1d) Table S1). Other details of the HBs are given in Table S1.
To study these effects quantitatively, pair wise RMSD values were computed by superposing each snapshot (coordinates saved at every 1 ns of the last 10ns of simulations) of the MD trajectory of DapB-N-P with the corresponding snapshots of DapB-Hyb. A similar matrix was also generated by superposing only the binding site residues of both the systems. These matrices (Fig 4a) indicated that there are significant structural differences between the two systems, especially in the binding site (C1) region. The L7 loop, which forms the C-terminal side of the binding site cavity C1 elevates about 3-4 Å in presence of the HM (Fig 4b) and interactions showed more than 10% occupancy (highlighted in Table S1, Fig 4c). The most stable interactions made by the ADP part of NADH were with D33 and K11 while PDC made stable interactions with H132, H133 and K136 of C1 in DapB-N-P.
It was found that the interactions formed by the nicotinamide part of NADH with A102, P103 and F105 are not stable (occupancy 0.06%, 0.14% and 0.60% respectively) throughout the simulation of DapB-N-P, although they are present in the crystal structure. In DapB-Hyb, out of 28 (Table S1, Figs S2 and S3), 11 interactions were found to be exhibiting more than 10% occupancy during the simulation (highlighted in Table S1, Fig 4d). The e-pharmacophore approach uses the energy components of the receptorligand binding from Glide XP 25 to generate structure-based pharmacophore models.
The models comprise six types of chemical features, viz., HB acceptor (A), HB donor  features. Considering the striking differences in the number and spatial orientations of the pharmacophoric features in the three types of models, we screened a huge chemical space using these models and compared the structures, druglike properties and target binding affinities of molecules screened by each type of models.

Comparative analysis of the molecules screened by the three types of pharmacophore models
As the three types of pharmacophore models were significantly different, it was interesting to study the similarities and diversities of molecules screened by each of them (Set-N, Set-P and Set-H compounds as described in the method). As NADH is a very abundant metabolite in human, the N-type pharmacophore models generated from NADH may not screen compounds with high specificity for Mtb-DapB. At the same time 2,6-PDC being a very small molecule, gives P-type pharmacophore models with only 3 to 4 features, which may fail to screen Mtb-DapB specific molecules. As observed from our MD simulation trajectories, interactions of the HMs with the hinge region residues bring about significant structural alterations in the binding site region, and none of the N-and P-types of model possesses features representing the hinge region interactions. So, the whole idea of generating the H-type models was to test them as tools to screen new compounds which not only bind to Mtb-DapB with high affinity and specificity, but also target the hinge region residues to alter the protein structure. The 133 P-set compounds were mostly single ring and low molecular weight fragment like structures. So, they were excluded in the comparative analysis of properties considering their small size and lack of chemotype diversity. The N-type models screened 982 molecules while the H-type models screened 387 compounds.
As a preliminary comparison, the structural similarity between Set-N/Set-P compounds were evaluated as Tanimoto coefficient (TC) [34] using OpenBabel [35] and a small inhouse shell script. The results are presented as a color-coded matrix (Fig 6a), where a TC value from 0 (lowest) to 1 (highest) is coded by red through yellow to green. Thus, the red shades represent lower similarity and the yellow to , the average score for the H-set compounds is slightly higher than the N-set compounds. As most of the scoring functions are biased towards larger ligands, we also considered the ligand efficiency indices [37] (Figs 6c and S6), which effectively relatively more solvent exposed groove like region (where the ADP part of NADH binds in the crystal structure), the H-set compounds occupy the C1 cavity (Fig 6b). this number is as low as 4% in case of N-set compounds (Fig 6c, 3 rd graph). This observation clearly indicates that the H-type pharmacophores are able to screen compounds that specifically interact with the hinge region residues in the C1 cavity and hence may impart better enzyme inhibitory effect.
Apart from the binding specificity, the druglikeliness of the N-and H-set compounds were also compared to further prioritize the H-type pharmacophore models. Fig S7 shows  compounds as compared to the H set compounds (Fig 6d) showing better druglikeliness of the later. The #metab descriptor is a predicted value representing number of likely metabolic reactions gives an estimation of the off-target interactions and toxicity of the compounds. The N set compounds show a higher number as compared to the H set compounds. Hence with all these comparative observations of the structural and physicochemical properties, and drug-likeliness scores, we can summarize that, the hybrid pharmacophore models lead to a structurally diverse and more druglike chemical space.

Application of the hybrid pharmacophore models to decipher the mechanism of action of anti-Mtb compounds
Considering the abilities of the hybrid pharmacophore models, they can have versatile applications. For example, they can be directly implemented to screen the huge chemical space to suggest potential binders of Mtb-DapB, which may act as start points for inhibitors design, the spatial locations of the pharmacophoric features may also exploited for fragment based de novo drug design and scaffold hopping.

Fig 7.
Interactions of two of the MA library compounds with the C1 cavity residues (please refer the legend given in Fig 5). The respective Glide XP docking scores are mentioned in the brackets.
In this study we demonstrate the application of the models to screen a library of known antitubercular compounds to find their possible mechanism of action. The MA library was created by extracting 2950 Mtb active compounds from ChEMBL and also manually from literature. These molecules have been reported to have ability to kill Mtb at concentration < 2 µM in the whole cell based essays, but their mechanism of action is not known. Screening this library with the eight H-type models returned a unique set of 31 compounds, which match five or more pharmacophoric features.
These compounds were docked to the DapB-Hyb snapshots to estimate their binding potentials as well as analyze their interactions with the crucial residues identified from the MD studies. Eight compounds were found to have the XP docking score above 7 (Table S4). These molecules make H-bonds, π -π, cation-π and salt bridge interactions with both the hinge region as well as other C1 residues (Fig 7, Table S4).

Hence, we predict these molecules might be acting on Mtb cells through inhibition of
Mtb-DapB by competitively binding to both the cofactor and substrate binding sites and/or restricting the inter-domain motions. However, these molecules need further experimental investigation to confirm their mechanism of action. These newly identified interactions may be utilized as Mtb-DapB specific screening filters and provide new insight for structure activity relationships. We report these interactions in the form of hybrid e-pharmacophore models, which can be used as key filters in structure based de novo design/screening of Mtb-DapB inhibitors. The compounds screened by the hybrid models were found to specifically make interactions with the hinge region residues and have reasonably good binding affinities with Mtb-DapB as estimated from the docking studies. Also, the hybrid pharmacophore models are capable of screening compounds which are highly diverse, structurally different and found to have better druglike properties than those screened by the models generated from NADH and 2,6-PDC individually. The study also demonstrates one of the many applications of these models by screening a library of compounds which are Mtb membrane permeable and show MIC values below 2µM to find if any of them acts through Mtb-DapB inhibition. Eight of them were found satisfying the pharmacophoric requirements of binding sites as well as making energetically favourable binding interactions. These findings may provide some clue for further experimental investigation to understand the mechanism of action of these compounds. Considering its broad applicability, we propose the hybrid epharmacophore modeling as a novel and effective strategy to venture into novel and potential anti-Mtb-DapB chemical space. Though in this article, we have discussed the dynamics-based hybrid pharmacophore strategy taking Mtb-DapB as an example, it can also be tested with a broad range of other drug targets.