Novel phytochemical–antibiotic conjugates as multitarget inhibitors of Pseudomononas aeruginosa GyrB/ParE and DHFR

Abstract

Background

There is a dearth of treatment options for community-acquired and nosocomial Pseudomonas infections due to several rapidly emerging multidrug resistant phenotypes, which show resistance even to combination therapy. As an alternative, developing selective promiscuous hybrid compounds for simultaneous modulation of multiple targets is highly appreciated because it is difficult for the pathogen to develop resistance when an inhibitor has activity against multiple targets.

Methods

In line with our previous work on phytochemical–antibiotic combination assays and knowledge-based methods, using a fragment combination approach we here report a novel drug design strategy of conjugating synergistic phytochemical–antibiotic combinations into a single hybrid entity for multi-inhibition of P. aeruginosa DNA gyrase subunit B (GyrB)/topoisomerase IV subunit B (ParE) and dihydrofolate reductase (DHFR) enzymes. The designed conjugates were evaluated for their multitarget specificity using various computational methods including docking and dynamic simulations, drug-likeness using molecular properties calculations, and pharmacophoric features by stereoelectronic property predictions.

Results

Evaluation of the designed hybrid compounds based on their physicochemical properties has indicated that they are promising drug candidates with drug-like pharmacotherapeutic profiles. In addition, the stereoelectronic properties such as HOMO (highest occupied molecular orbital), LUMO (lowest unoccupied molecular orbital), and MEP (molecular electrostatic potential) maps calculated by quantum chemical methods gave a good correlation with the common pharmacophoric features required for multitarget inhibition. Furthermore, docking and dynamics simulations revealed that the designed compounds have favorable binding affinity and stability in both the ATP-binding sites of GyrB/ParE and the folate-binding site of DHFR, by forming strong hydrogen bonds and hydrophobic interactions with key active site residues.

Conclusion

This new design concept of hybrid “phyto-drug” scaffolds, and their simultaneous perturbation of well-established antibacterial targets from two unrelated pathways, appears to be very promising and could serve as a prospective lead in multitarget drug discovery.

Keywords:

• hybrid compounds
• multi-target inhibition
• drug resistance
• dihydrofolate reductase
• DNA gyrase subunit B
• topoisomerase IV subunit B
• rational drug design

Introduction

Pseudomonas aeruginosa is one of the most opportunistic and challenging pathogenic bacteria, with constant evolution of resistance and the complexity of multidrug resistant phenotypes due to the extensive use of antibacterial agents by monotherapeutic strategy.¹ Multidrug resistant phenotypes in P. aeruginosa occur through the acquisition of multiple imported resistance mechanisms, coupled with chromosomally encoded resistance mechanisms, accumulation of multiple chromosomal changes over time, and/or a single mutational event leading to overexpression of one or more efflux pumps.² This continuing trend of developing drug resistance can severely limit the therapeutic options for treatment of serious infections caused by pathogens, specifically, P. aeruginosa. Alarmingly, the prospects for novel antipseudomonal drugs in the next few years are quite poor, with most of the drugs in the pipeline being more or less new derivatives of existing families of antimicrobial agents.³ As an alternative strategy, antibacterial combination therapy is being used as a standard practice in many hospitals for a broader empirical coverage by simultaneously administering two different antibacterial agents, leading to a synergistic effect. However there are no clinical data supporting its effectiveness against gram negative infections as well as the issues of drug toxicity, side effects, and selection of multidrug resistant organisms.^4,5 Taken together, it is imperative to identify novel chemical entities through effective therapeutic strategy to combat the threat of multidrug resistance (MDR).

Antibacterial polypharmacology is an emerging strategy that overcomes the limitations of both monotherapy and combination therapy by targeting multiple proteins simultaneously in the disease-associated network.^6–8 Recently, there has been greater recognition that multi-action hybrid compounds, through fragment combination, can be a potent strategy against robust pathogens because of their ability to delay the onset, decrease the likelihood, and overcome the existing resistance mechanisms with reduced toxicity concerns.⁶ In a fragment combination approach, two dissimilar entities are combined together with a noncleavable linking group designed to bind independently to two different biological targets and synchronously accumulate at both the target sites.^9,10 Furthermore, in comparison to a combination of drugs, multitarget inhibitors have greater predictable pharmacokinetic and pharmacodynamic relationships, as they involve administration of a single drug.¹¹ This fragment combination strategy holds significant promise and has been successfully applied in the design and development of novel, dual-action, hybrid antibacterial agents, such as hybrids, that delay the emergence of bacterial resistance by inhibiting bacterial protein synthesis and fluoroquinolone targets,¹² berberine (plant alkaloid, sensitive to NorA efflux) – (2-phenyl-5-nitro-1H-indole, NorA efflux pump inhibitor) hybrid antibacterial, insensitive to MDR efflux,¹³ trimethoprim–fluoroquinolone hybrids for dual action in dihydrofolate reductase enzyme and fluoroquinolone targets (DNA gyrase subunit A, [GyrA] and topoisomerase IV subunit A, [ParC]).¹⁴ Thus, the evidence increasingly shows that discovery of druglike entities with predefined promiscuity is emerging as a productive area of research.

In line with this, we employed a combined knowledge-based and fragment combination approach, focused on the design and development of novel multitarget hybrid compounds that can act by targeting GyrB/ParE and DHFR enzymes through optimization of phytochemicals by conjugating with antibiotics. To the best of our knowledge, this study is the first to successfully combine knowledge from the scientific literature with previous drug combination data and a fragment combination approach towards the identification of novel multi-targeting “phyto-drug” conjugates. In silico-based approaches have always shown promising potential in facilitating drug discovery, and these methods have also been successfully employed for the search and design of selective multitarget compounds.¹⁵ Furthermore, we have integrated various effective computational tools to investigate the molecular, electronic, and pharmacokinetic properties, binding affinity and stability of the inhibitor–enzyme complexes of our designed compounds.

Methods and computational details

Ligand data set preparation and properties predictions

Marvin sketch (ChemAxon, Budapest, Hungary) was used to create the SMILES (.smi) format for the designed inhibitors. Open Babel online interface (http://openbabel.sourceforge.net/), through shell script and using .smi as input, was used to generate primary 3D structures, which were then saved in SDF file format.¹⁶ The structures were subjected to conformational search calculations in gas phase using the Merck molecular force field method in the SPARTAN’10 program (WaveFunction, Inc, Irvine, CA, USA). The minimum energy conformation was chosen for geometry-optimization using the Austin model 1 (AM1) semiempirical algorithm. The optimized structures were reoptimized using the density functional theory hybrid, B3LYP/6-31G** basis set functional.¹⁷ The stereoelectronic properties such as HOMO (highest occupied molecular orbital) and LUMO (lowest unoccupied molecular orbital) energies, band gap (HOMO–LUMO), MEP, dipole moment and polarizability, physicochemical properties (molecular volume, polar surface area [PSA]), and absolute energies were calculated from single-point energy calculations on the complete geometry-optimized structures of the designed compounds using B3LYP/6-31G** basis set. In addition, percentage of absorption (%ABS) was calculated by using the following equation,

$\%\text{ABS }=109-(0.354\times \text{ PSA})$(1)

based on the computed PSA values.¹⁸

To obtain additional insight on the toxicity risks (mutagenicity, tumorigenicity, irritation, and reproductive effects) and physicochemical properties (molecular weight, cLogP, solubility, PSA, number of hydrogen bond donors and acceptors, number of rotatable bonds, drug-likeness, and overall drug score) of the designed hybrid compounds, GyrB/ParE inhibitors (cyclothialidine and novobiocin) and DHFR inhibitors (methotrexate and trimethoprim), phytochemicals (protocatechuic acid and gallic acid), and sulfonamides (sulfamethoxazole and sulfadiazine), were calculated using OSIRIS property explorer¹⁹ and Molinspiration property prediction toolkit (Molinspiration Cheminformatics, Slovensky Grob, Slovak Republic) (http://www.molinspiration.com/cgi-bin/properties). The toxicity risk assessment indicated potential risks in the drawn structure concerning the risk category specified, and the process relied on the precomputed set of structural fragments that give toxicity alerts if they were encountered in the drawn 2D structures. Properties with high risks of undesired effects, like mutagenicity or a poor intestinal absorption, were shown as a negative sign, whereas a positive sign indicated drug-conform behavior.

Molecular docking

The MolDock docking engine in Molegro virtual docker (CLC bio, Aarhus N, Denmark) was used for studying the binding modes of the designed hybrid compounds in the active site cavities of the modeled P. aeruginosa GyrB/ParE and DHFR enzymes (see Supplementary materials), which was based on a new heuristic search algorithm (MolDock score) that combines differential evolution with a cavity prediction algorithm.²⁰ In our docking experiments, a MolDock grid scoring function using template docking with default values: −500 overall strength and 0.4 Å energy grid resolution was used to evaluate the energy between the ligand and the target enzyme. Grid resolution, number of runs, population size, maximum iterations, pose generation energy threshold, simplex evolution max steps, and neighbor distance factor were set as 0.30 Å, 20, 50, 1500, 100, 1.00 for each run, respectively, using the MolDock SE algorithm. The ligands from the crystal structures of Escherichia coli GyrB/ParE and DHFR were transferred into the workspace, keeping the orientation as a control and were kept as the reference ligand. The complete geometry-optimized structures of the hybrid compounds and the generated protein homologs were also transferred, and hydrogen molecules were added to both ligands and protein molecules using the preparation wizard in the Molegro workspace. During import of the 3D structures of the ligands, charges and bond orders were assigned, the torsional angle of the 3D structures was also determined, and all acyclic single bonds were set as flexible. Binding sites in the electrostatic surface of the protein were identified using the grid-based cavity prediction algorithm. A total of five cavities were detected, the prepositioned reference ligand in the active site cavity was identified, and the docking was constrained to the predicted active site cavity. Multiple poses were returned for each run with the root mean square deviation (RMSD) threshold set to 1.00 Å. The pose with the highest rerank and MolDock score was retained in the workspace for detailed evaluation of the ligand binding at the active site cavity. The rerank score uses a weighted combination of the terms used by the MolDock score mixed with a few additional terms (the rerank score includes the steric terms which are Lennard–Jones approximations to the steric energy; the MolDock score uses a piece-wise linear potential to approximate the steric energy).²⁰ The rerank scoring function improved the docking accuracy by identifying the most promising docking solution from the solutions obtained by the MolDock docking algorithm.²⁰ The rerank score provided an estimate of the strength of the interaction. It was not calibrated in chemical units, and it did not take complex contributions such as entropy into account. Even though the rerank score might be successful in ranking different poses of the same ligand, it might be less successful in ranking poses of different ligands.²⁰ Along with both MolDock and reranking scores, we also predicted binding affinities using a calibrated model that is included in the Molegro virtual docker. The binding affinity model was trained using a data set of more than 200 structurally diverse complexes from Protein Data Bank (PDB) with known binding affinities.²¹ Hence, in our docking experiments we used this recommended strategy of ranking the docking results by their rerank scores and subsequently the binding affinity measure to get high ranked poses. The validation of the docking protocol was carried out by redocking the imported reference ligands from their respective experimental PDB structures in the predicted active site cavity of the model using the RMSD measure. The RMSD (Å) values for the redocked inhibitor methotrexate (MTX) and substrate NADP in the modeled P. aeruginosa DHFR active site were 0.3968 Å and 0.5794 Å, respectively, wherein the substrate, ADPNP (5′–adenylyl-β-γ-imidodiphosphate, the nonhydrolyzable analog of ATP) in the P. aeruginosa GyrB/ParE were 0.2529 Å/0.2052 Å, indicating high similarity between predicted and experimental binding mode. Furthermore, the binding conformations of ligands in the enzyme active site cavity were analyzed by visualizing the H-bond and electrostatic interactions formed between the ligands and the active site residues.

Molecular dynamics simulations

The molecular dynamics simulations for the corresponding substrate bound enzymes and top ranked docked inhibitor in both the enzymes were performed with the YASARA dynamics package version 10 (YASARA Biosciences, Vienna, Austria).²² A periodic simulation cell of at least 20 Å larger than the protein was employed with explicit solvent around the P. aeruginosa enzyme complexes (I–GyrB, IV–DHFR VII–ParE, II–ADPNP–GyrB, V–MTX-NADP-DHFR, VIII–ADPNP-ParE), and top ranked hybrid compound-bound enzyme complexes (III-phyto-drug [PD]_2a-GyrB, VI–PD_2a-NADP-DHFR, and IX–PD_2a-ParE) were energy minimized to correct the covalent geometry and to remove bumps. The AMBER03 force field was used with Van der Waals pairs cut-off distance 7.86 Å and long range statistics algorithm calculated using the Particle Mesh Ewald method.²³ The automatic force field parameter assignment was carried out by AutoSMILES,²⁴ which assigned pH-dependent fractional bond orders and protonation patterns. This was followed by geometry-optimization of the structures with the COSMO solvation model²⁴ using semiempirical AM1 Mulliken point charge calculations and assignment of AM1BCC (AM1 bond charge correction),²⁵ atom, and bond types, with further refinement using known restrained electrostatic potential charges. And finally, GAFF (General AMBER force field)²⁶ atom types and remaining force field parameters were assigned. The hydrogen bonding network of the complexes were optimized by a method that positions polar hydrogen atoms in protein structures by optimizing the total hydrogen bond energy.²⁷ The simulation cell was filled with TIP3 water, and Na⁺ and Cl⁻ counter ions were added at the most favorable positions to neutralize the ion strength in the cell. A few water molecules were deleted to readjust the solvent density to 0.997 g/mL, and the pressure was controlled by rescaling the simulation box dimensions to maintain the water density. The pKa values were predicted using the empirical pKa prediction equation that is approximated as a function of electrostatic potential, hydrogen bonds, and accessible surface area by Ewald summation implemented in YASARA.²⁸ A short molecular dynamics simulation was run on the solvent only, and the entire system was energy minimized by steepest descent minimization to remove conformational stress. This was followed by simulated annealing minimization using a time step of 2 femtoseconds and atom velocities scaled down by 0.9 every 10th step until convergence (energy improved by less than 0.1% during 200 steps). The molecular dynamics simulation was then initiated, with temperature fixed at 298 K and multiple time steps set for intramolecular and intermolecular forces at 1.33 femtoseconds and 4 femtoseconds, respectively. Each complex was subjected to 3000 picoseconds (ps) of molecular dynamics simulations, and snapshots were saved every 1 ps for data analysis. The trajectories of the molecular dynamics simulations were analyzed for the equilibrium stability by measuring the RMSD of the complexes and the root mean square fluctuation (RMSF) of residues around the ligand-binding active sites.

Results

Multi-action hybrid compounds design hypothesis

As a first step, we surveyed the biomedical literature to establish a link between the targets (GyrB/ParE and DHFR) based on the fact that if the same compound has activity against different targets they might be related (Figure 1B).²⁹ It has been reported that the plant flavonoid, epigallocatechin gallate and its derivatives access the active site cavities of GyrB, ParE and DHFR) enzymes,^30–32 which is illustrated in Figure 1B. In support of this, we proposed earlier that phytochemicals (structural analogs of epigallocatechin gallate, protocatechuic acid and gallic acid) are synergistic in combination with the antifolate drugs, sulfamethoxazole and sulfadiazine (which target dihydropteroate synthase enzyme in the folic acid biosynthesis pathway) by binding to DHFR, and are additive in combination with ciprofloxacin (which targets GyrA and ParC) by binding to GyrB and ParE against P. aeruginosa isolates^33,34 (Figure 1A). It has also been reported that derivatives of benzene sulfonamide groups also bind to the target enzymes,^35,36 which is illustrated in Figure 1B. Additionally, important interactions for the inhibition of the enzymes GyrB/ParE and DHFR of E. coli have been identified based on their individual pharmacophore models from molecular dynamics simulations.^37,38 Using these pharmacophore models we derived a common pharmacophoric feature model for multiinhibition of these enzymes (Figure 1B), which are (1) an aromatic element with hydrogen-bond donor elements mapping the adenine moiety of ATP in GyrB/ParE and the pterin/pteridine moiety of dihydrofolate/methotrexate in DHFR,^37,38 and (2) ribose and phosphate groups of ATP in GyrB as well as p-aminobenzoate linker and glutamate tail region of folate, mapping for aromatic/hydrophobic moiety with hydrogen-bond acceptor elements.^37,38 Thus, we embarked on using the simple phenols (gallic acid and protocatechuic acid), which are hydroxyl-substituted aromatic ring compounds identified from our previous work, to satisfy the common pharmacophoric feature: (1) conjugated with a noncleavable linking group (hydrophobic element) to another aromatic element with hydrogen bond acceptor moieties and (2), the sulfonamide groups (sulfamethoxazole and sulfadiazine) for the design of multitargeting “phyto-drug” hybrid compounds (Figure 1C).

Figure 1 How did we connect the dots? Hypothesis and design of multitargeting compounds. (A) Drug–phytochemical interactions identified using combination assays. Green arrows signify synergistic interaction of phytochemicals (protocatechuic acid and gallic acid, represented as blue arrowhead) in combination with sulfonamide antibiotics (sulfamethoxazole and sulfadiazine, represented as violet arrowhead). Sulfonamides block the dihydropteroate synthase enzyme, while the phytochemicals target the subsequent enzyme dihydrofolate reductase (DHFR) in the folate biosynthesis pathway leading to a synergistic effect. Red arrows signify indifferent interaction of phytochemicals (protocatechuic acid and gallic acid, represented as blue arrowhead) in combination with fluoroquinolone antibiotics (ciprofloxacin, represented as orange arrowhead). Ciprofloxacin binds to DNA gyrase subunit A (GyrA) and topoisomerase IV subunit A (ParC), while phytochemicals target the subunit B (GyrB and ParE) of both these enzymes leading to an indifferent effect via equivalent action of DNA replication inhibition; (B) Knowledge-based approach was used to establish relationship between therapeutic targets (GyrB/ParE and DHFR) identified based on our previous drug–phytochemical interaction network. Mechanism of action of epigallocatechin gallate (plant flavonoid, red arrowhead) and previously reported sulfonamide derivatives targeting GyrB/ParE and DHFR enzymes, gives an illustration of the ability of these compounds to bind to the identified target sets. Derived common pharmacophoric features for multi-inhibition of GyrB/ParE and DHFR enzymes based on the previously reported individual pharmacophore models; (C) Schematic illustration of the fragment combination method used to derive the hybrid compounds, phytochemicals (blue arrowhead) linked with sulfonamide groups (violet arrowhead) using noncleavable linker (orange line) that can simultaneously bind to GyrB/ParE and DHFR enzymes.

In silico calculations of physicochemical, stereoelectronic properties, drug likeness, and toxicity risks

In order to identify the ultimate drug candidate it is first necessary to start with right lead compounds of drug likeness and ligand efficiency.³⁹ Thus in recent years, various in silico-based predictions of physicochemical properties with the consideration of drug likeness have gained importance due to their potential to rapidly eliminate drug candidates prior to synthesis. It is also evident that several physicochemical properties are useful indicators for drug likeness of a compound with an impact on specific pharmacokinetic processes like permeation, intestinal absorption, reversible binding, distribution, and excretion of a compound.⁴⁰

As a preliminary step in our study, we investigated relevant toxicity risks, physicochemical properties, bioavailability, drug likeness, and overall drug scores of the designed hybrid compounds using various in silico-based prediction toolkits. The predicted toxicity risks shown in Table 1 indicate that all the hybrid compounds are nonmutagenic, nontumorigenic, nonirritating, nonreproductive, and are comparable to standard traded drugs. For comparison, we evaluated the toxicity risks for commercially available GyrB/ParE inhibitors (cyclothialidine and novobiocin) and DHFR inhibitors (methotrexate and trimethoprim) of which cyclothialidine and trimethoprim showed nontoxic effects, whereas novobiocin and methotrexate showed risks concerning reproductive and tumorigenic effects, respectively. Also, the predicted toxicity risks for phytochemicals (protocatechuic acid and gallic acid) showed mutagenic effects for protocatechuic acid, while gallic acid showed both mutagenic and reproductive effects. Among sulfonamides (sulfamethoxazole and sulfadiazine), sulfamethoxazole showed tumorigenic risk effects, and sulfadiazine was found to be associated with a reproductive toxicity alert. Overall, the proposed hybrid structures (PD_1a, PD_1b, PD_2a, and PD_2b) in Figure 2 showed nontoxic behavior with highly desirable physicochemical parameters disclosing their potential as promising agents for antimicrobial therapy.

Table 1 Calculated toxicity risks, molecular properties, drug likeness, and overall drug score of the compounds investigated

Compounds	Toxicity risks prediction				Molecular properties
	MUT	TUM	IRR	REP	MW	cLP	S	nON	nOHNH	nViolations	nRotb	D_L	D_S
PD_1a					449	1.29	−3.67	11	5	1	9	4.2	0.72
PD_1b					446	0.99	−2.95	11	5	1	9	4.81	0.77
PD_2a					433	1.59	−3.97	10	4	0	9	3.75	0.71
PD_2b					430	1.29	−3.24	10	4	0	9	4.33	0.77
Protocatechuic acid					154	0.9	−1.04	4	3	0	1	−0.12	0.43
Gallic acid					170	0.6	−0.74	5	4	0	1	0.12	0.27
Sulfamethoxazole					253	0.7	−3.02	6	3	0	3	3.31	0.54
Sulfadiazine					250	0.41	−2.29	6	3	0	3	4.09	0.56
Cyclothialidine					625	−2.62	−2.31	16	9	3	6	2.14	0.54
Novobiocin					598	2.81	−4.93	13	6	3	9	0.97	0.23
Methotrexate					454	−0.69	−3.77	13	7	2	9	−7.09	0.22
Trimethoprim					290	1.23	−3.32	7	4	0	5	4.95	0.87

Notes: , Toxic free; , Toxic behavior.

Abbreviations: D_L, druglikeness; D_S, drug score; IRR, irritating effects; cLP, cLogP; MUT, mutagenicity; MW, molecular weight; nON, number of hydrogen bond acceptors; nOHNH, number of hydrogen bond donors (OH and NH groups); nRotb, number of rotatable bonds; REP, reproductive effects; S, Solubility; TUM, tumorigenicity.

Figure 2 Structures of the designed hybrid compounds, “phyto-drug” (phytochemicals–antibiotics conjugates, (A) PD_1a; (B) PD_1b; (C) PD_2a; (D) PD_2b.

The Lipinski’s RO5, used as a filter for drug likeness, defines four simple physicochemical parameter ranges (molecular weight ≤ 500, cLogP ≤ 5, H-bond donors ≤ 5, H-bond acceptors ≤ 10) associated with acceptable aqueous solubility and intestinal permeability. Molecules violating more than one of these rules may have problems with oral bioavailability.⁴¹ The physicochemical properties calculated, such as molecular weight of the proposed compounds are below 450 (Table 1), and thus are more likely to have higher absorption. The cLogP value is the logarithm of its partition coefficient between n-octanol and water, measuring the compound’s hydrophilicity¹⁹ The estimated cLogP values for the proposed compounds are less than 5.0, indicating high lipophilicity and higher absorption and permeation. The aqueous solubility (log S) measure significantly affects the drug absorption/distribution characteristics, and is a unit stripped logarithm of the solubility (moles/L). It is suggested that more than 80% of the drugs in the market have aqueous solubility values greater than −4.¹⁹ The solubility measures for the proposed compounds in Table 1 are in the expected range of −2.95 to −3.97 and are comparable to that of the standard commercially available drugs. The number of Lipinski’s RO5 violations is one each for designed hybrid compounds PD_1a, PD_1b and none for PD_2a and PD_2b, which are within the range sets defined for oral bioavailability. The number of rotatable bonds and PSA are now widely used filters for drug likeness, and those compounds which meet only the two criteria of: (1) ≤ 10 rotatable bond counts, and (2) PSA ≤ 140 Ų (or ≤ 12 H-bond donors and acceptors) will have increased oral bioavailability in the rat.⁴² The calculated rotatable bond counts (Table 1) of the proposed compounds are equal to 9, and PSA are within the range for compounds PD_2a and PD_2b, and just above for compounds PD_1a and PD_1b (Table 2). On the other hand using PSA, we calculated the %ABS according to equation (1)$\%\text{ABS }=109-(0.354\times \text{ PSA})$(1) : %ABS = 109 − (0.354 × PSA) as reported by Zhao et al,¹⁸ the calculated absorption percentages are in the range of 55%–65% for all the proposed compounds (Table 2). Drug likeness scores are calculated by summing up score values, partially based on topological descriptors, and fingerprints of Molfile MDL (Accelrys, San Diego, CA, USA) structure keys including properties, such as cLogP and molecular weight, are estimated to be in the positive range for 80% of the traded drugs.¹⁹ The drug likeness scores for the proposed compounds are positive (Table 1), indicating that the molecules predominantly contain fragments which are frequently present in commercial drugs, and hence have the possibility of being ‘drug-like’. Finally, the overall drug score (Table 1), combining drug likeness, cLogP, aqueous solubility, molecular weight, and toxicity risks describe the potential to qualify for a drug on 0 to 1 scale.¹⁹ The overall drug scores were in the range of 0.71–0.77 (Table 1), indicating the proposed compounds potential to qualify as suitable drug candidates.

Table 2 Calculated absolute energies, stereoelectronic, and molecular properties at B3LYP/6-3IG** basis set levels, and percentage of absorption of the designed hybrid compounds

Compounds	Energy (in hartress)	Energy (aq) (in hartress)	Energy (S) (KJ/mol)	%ABS	Volume (Å^Citation3)	PSA (Å^Citation2)	EHOMO (eV)	ELUMO (eV)	Band gap (ΔE) (eV)	Dipole moment (debye)	Polarizability
PD_1a	−1899.6	−1899.66	−109.8	55	402.14	152.39	−6.26	−1.38	−4.88	7.35	72.84
PD_1b	−1878.6	−1878.64	−127.37	59	400.06	142.42	−5.84	−1.49	−4.35	7.71	72.8
PD_2a	−1824.4	−1824.44	−94.78	65	394.08	125.09	−6.16	−1.13	−5.03	7.06	72.15
PD_2b	−1803.3	−1803.42	−124.07	63	393.39	129.05	−6.05	−1.25	−4.8	5.21	72.15

Abbreviations: %ABS, absorption; ΔE, ΔE_HOMO-LUMO; Energy, total energy in vacuo; Energy(aq), total energy in aqueous medium; Energy(S), aqueous salvation energy; PSA, polar surface area.

In order to gain more insight into the molecular recognition through ligand-receptor interactions and structural factors that are involved in the multispecific activity of the designed hybrid compounds, stereoelectronic properties such as molecular orbital energies (HOMO and LUMO), band gap (difference between HOMO and LUMO), and MEP profiles were calculated using the B3LYP/631G** density functional theory on the geometry-optimized structures. Calculated ab initio quantum chemical B3LYP/6-31G** calculations and electronic properties of the compounds are listed in Table 2. HOMO and LUMO orbital energies are related to ionization potential and electron affinity and their respective frontier orbitals are associated to the molecule’s reactivity ie, HOMO energy is susceptible to electrophilic attack (donate e⁻) and LUMO energy susceptible to nucleophilic attack (accept e⁻).⁴⁴ The reactivity index (band gap) of the compounds with small difference implies high reactivity, and large difference implies low reactivity in reactions.⁴⁵ The energy gap (HOMO – LUMO, shown in Table 2) of the proposed compounds was found to be within a narrow range of −4.35 to −5.03, indicating that these energies permit electron transfer and exchange making those compounds very highly reactive and similar in nature.

Three dimensional HOMO and LUMO isosurfaces over the Van der Waals surface and the MEP mapped onto total electron density surface of all the compounds are shown in Figure 3A–C. For the three dimensional HOMO orbitals of all the compounds (Figure 3A), the predominant HOMO (electron-rich) orbitals were consistently observed on the hydroxyl substituted aromatic ring of the phytochemicals (gallic acid and protocatechuic acid), indicating strong electrophilic affinity, and the most significant LUMO (electron-poor) orbitals (Figure 3B) of all the compounds were consistently observed on the sulfonamide rings, indicating the localization of strong electron acceptor elements, and is thereby susceptible to strong nucleophilic attack. Additionally, the MEP is a useful descriptor for understanding electrophilic/nucleophilic sites as well as hydrogen bonding interactions. In general, the negative potential sites (red color in Figure 3) of MEP represent regions of electrophilic reactivity and interaction through π-π bonding with the aromatic systems of interacting enzyme residues, while positive potential sites (blue color in Figure 3) represent regions of nucleophilic reactivity.⁴⁶ As expected, the MEP maps of all the compounds (Figure 3C), clearly show the presence of both a positive potential site and a negative potential site around the hydroxyl groups in the aromatic ring of the phytochemicals, suggesting that these hydroxyl groups can act as both hydrogen bond donor and acceptor elements, indicating their participation in hydrogen bond interactions with hydroxyl and carbonyl groups of key active site residues in GyrB (Asp75, Thr167, Tyr111, and water molecule), ParE (Asp68, Glu45, and water molecule) and DHFR (Asp30, Ile8, Tyr110, and Ile104) enzymes. The predominant red regions (negative potential sites) in the MEP profiles of all the compounds are localized around the sulfonyl functional groups (oxygen atoms) of the sulfonamide rings, nitrogen (N7) and oxygen (O8) atom in the oxazole ring of sulfamethoxazole, nitrogen (N7) atoms in the pyrimidine ring of sulfadiazine, oxygen (O8) atom and the carbonyl oxygen (O8) atom in phytochemicals, suggesting that these oxygen and nitrogen atoms can act as hydrogen bond acceptors and possibly participate in hydrogen bonding interaction with donor groups of key active site residues in GyrB (Tyr111, Asn48, Val122, Lys105, Leu117, His118, and Gly119), ParE (Tyr104, Asn41, Ile115, Gly114, Gly112, His111, and Leu110), and DHFR (Arg55 and Arg60) enzymes. Also, regions at the aromatic rings in all the compounds show red regions with decreasing intensity, indicating weak positive electrostatic potential and thereby implying a hydrophobic nature of the aromatic rings. The above mentioned molecular orbitals and MEP profile features appear to be consistent with the key pharmacophoric derivatives required for potent multi-inhibition of GyrB/ParE and DHFR enzymes.

Figure 3 Representations of (A) HOMO (B) LUMO and (C) Molecular electrostatic potential (MEP) isodensity surfaces calculated at B3LYP/6-31G** basis set levels on the optimized geometry of the hybrid compounds.

Molecular docking

To better understand the binding modes of the designed hybrid compounds in both GyrB/ParE and DHFR enzymes, molecular docking studies were carried out using the Molegro virtual docker program. The docking results of all the proposed compounds along with ADPNP (GyrB/ParE substrate)/methotrexate (DHFR inhibitor) in P. aeruginosa GyrB/ParE and DHFR are ranked based on their binding affinity scores (Tables 3 and 4, and Table S1).

Table 3 Docking results of ADPNP and designed hybrid compounds (ranked on the basis of binding affinity scores) at the ATP binding site in P. aeruginosa GyrB

Compounds	Scoring function		EInter	EIntra	Hbond	LE1	LE3	Binding affinity (kJ/mol)	Residues involved in hydrogen bonding
	MolDock	Rerank
PD_2a	−202.69	−164.34	−208.33	5.63	−23.64	−6.76	−5.48	−36.23	Asn48, Val122, Tyr111, Ser49, Asp75, Gly79, Glu52, Thr167, Leu117, His118, Gly119, HOH
PD_1a	−202.19	−164.58	−209.77	7.58	−26.33	−6.52	−5.31	−35.23	Asn48, Val122, Tyr111, Gly79, Thr167, Ser49, Asp75, Lys105, Lue117, His118, Gly119, HOH
PD_2b	−178.49	−148.53	−204.5	26.01	−20.36	−5.95	−4.95	−34.37	Ser49, Asp75, Tyr111, Gly79, Glu52, Asn48, Val122, His118, Gly119, Val120, Gly121, HOH
PD_1b	−192.07	−154.37	−209.47	17.41	−24.76	−6.19	−4.98	−34.33	Ser49, Asp75, Thr167, Gly79, Lys105, Asn48, Val122, His118, Gly119, Val120, Gly121, Tyr111, HOH
ADPNP	−234.13	−183.03	−219.98	−14.15	−28.94	−7.55	−5.9	−33.13	Val120, Gly121, Lys339, Leu117, His118, Lys105, Gly119, Asn48, Val122, Gly104, Tyr111, Asp75, Thr167, HOH

Abbreviations: MolDock, docking score function by MolDock; Rerank, reranking score function; E_Inter, ligand–protein (GyrB) interaction energy; E_Intra, internal energy of the ligand; HBond, hydrogen bonding energy; LE1, Ligand efficiency 1 – MolDock score divided by Heavy atoms count; LE3, Ligand efficiency 3 – Rerank score divided by Heavy atoms count; ADPNP, (nonhydrolyzable analog of ATP) is GyrB substrate.

Table 4 Docking results of methotrexate and designed hybrid compounds (ranked on the basis of binding affinity scores) at the folate binding site in P. aeruginosa DHFR

Compounds	Scoring function		EInter	EIntra	Hbond	LE1	LE3	Binding affinity (kJ/mol)	Residues involved in hydrogen bonding
	MolDock	Rerank
PD_2a	−139.98	−113.99	−154.61	14.63	−20.98	−4.67	−3.79	−34.81	Arg55, Thr123, Asp30, Ile8, Ala9, Tyr110, Ile104, Arg60
PD_1a	−150.91	−118.01	−161.61	10.70	−25.50	−4.87	−3.81	−34.77	Arg55, Ile8, Ala9, Ile104, Tyr110, Asp30, Thr123, Arg60
PD_2b	−136.41	−106.68	−155.81	19.39	−16.96	−4.54	−3.56	−30.78	Asp30, Thr123, Arg55, Ile8, Ala9, Ile104, Tyr110
PD_1b	−126.06	−87.87	−151.43	25.37	−16.77	−4.07	−2.83	−30.47	Ala9, Asp30, Thr123, Ile8, Ile104, Tyr110, Arg55, Arg60
MTX	−151.34	−119.61	−174.19	22.85	−18.08	−4.48	−3.54	−33.05	Asp30, Ile8, Ile104, Tyr110, Arg55, Arg60

Abbreviations: MolDock, docking score function by MolDock; Rerank, reranking score function; E_Inter, ligand – protein (DHFR) interaction energy; E_Intra, internal energy of the ligand; HBond, hydrogen bonding energy; LE1, Ligand efficiency 1 – MolDock score divided by Heavy atoms count; LE3, Ligand efficiency 3 – Rerank score divided by Heavy atoms count; MTX, (methotrexate) is DHFR inhibitor.

Docking of ADPNP and proposed hybrid compounds at the ATP binding site of GyrB

The 43-kDa domain of P. aeruginosa GyrB has two distinct subdomains, the N-terminal subdomain comprising an eight-stranded β sheet and four α helices containing the ATP-binding site, and the C-terminal subdomain consisting of a four-stranded β sheet and four α helices.⁴⁷ The GyrB in complex with the ATP analog, ADPNP, shows that the conformation is stabilized by H-bonds with residues Asn48, Asp75, Gly104, Lys105, Tyr111, Leu117, His118, Gly119, Val120, Gly121, Val122, Thr167, Lys 339, and Ile80, Ile96, and Val122 forming favorable hydrophobic interactions with the adenine and ribose moieties (Figure 4A), which are similar to the crystal structure of the E. coli GyrB43-ADPNP complex. In addition to these interactions, another important H-bond network was formed by the residues Asp75, Gly79, and Thr167 and the adenine ring of ADPNP with a conserved water molecule.⁴⁷ The residues that interact with the adenine ring of the ATP include Asp75, Tyr111, and Thr167. The ribose and phosphate groups are stabilized by Asn48, Lys105, Leu117, His118, Gly119, and Asn48 is one of the key residues required for binding of inhibitors.³⁸ Hence, a common feature can be rationalized in that the majority of the inhibitors of GyrB maintain a hydrogen bond network with Asp75, a highly conserved water molecule, and Asn48 residues, which is very important for binding of any new inhibitor class.^38,47 The docking results show that the proposed hybrid compounds have good binding affinities at the ATP binding site of GyrB, maintaining contacts with all the key residues required for inhibition. The major residues participating in hydrogen bond interactions with the designed inhibitors in the GyrB active site are tabulated in Table 3. Considering the top ranked hybrid compounds PD_2a and 1a, both having similar binding conformations and interactions (Figure 5A), the hydroxyl moiety in the carboxylic group of the protocatechuic acid forms hydrogen bonds with Asp75 and Ser49 and the carbonyl moiety forms a hydrogen bond with Asn48, similar to the adenine ring interaction of ATP. The 3′-OH and 4′-O group of the protocatechuic acid in PD_2a forms additional hydrogen bond contacts with Thr167, Gly79, a conserved water molecule, and Tyr 111. The additional hydroxyl group in the gallic acid moiety of PD_1a forms hydrogen bond interactions with Tyr111 and Lys105 (Figure 6A). The phenyl ring of sulfamethoxazole in both PD_2a and 1a sits in the ribose binding region forming hydrophobic contacts with residues Ile80, Ile96, and Val122, while the sulfonyl group positioning similar to the α-phosphate group of ADPNP forms hydrogen bond acceptor contacts with Asn48 and Val122 residues, and the oxazole ring orienting in between the β and γ-phosphate groups of ADPNP forms similar hydrogen bond acceptor interactions with Leu117, His118, and Gly119 (Figures 5A and 6A). Compounds PD_2b and 1b were found to bind, preferably maintaining the important interactions in a similar way to PD_2a and 1a compounds (Figures 5B and 6B). However, in comparison to PD_2a and 1a, the interactions established by compounds PD_2b and 1b are weaker, and the main reason for this effect may be attributed to the pyrimidine ring positioned directly over the γ-phosphate moiety, thereby limiting the hydrogen bond interactions with Gly119 and His118 residues. Supporting our claim, it has previously been demonstrated that phenolic groups have favorable binding energy at the adenine binding pocket if ATP and in the phosphate group binding site, oxazole has favorable clustering compared to benzene groups.³⁸

Figure 4 Binding modes of ADPNP and methotrexate in P. aeruginosa GyrB and DHFR active site cavities obtained from docking studies. (A) The nonhydrolyzable analog of ATP (ADPNP) (yellow stick model) docked at ATP binding site of GyrB; (B) DHFR inhibitor methotrexate (yellow stick model) docked at folate binding site of DHFR and the cofactor NADP is shown as aqua stick model.

Notes: The residues are shown as wireframe colored by their element, and the hydrogen bonds are illustrated as dotted green lines. The water molecule in (A) is shown as small red colored ball.

Figure 5 Docking of designed hybrid compounds at ATP binding site of P. aeruginosa GyrB. (A) Binding pose of the top scored compound PD_2a (grey stick model); (B) PD_2a and PD_2b (stick model colored by their element) superimposed with ADPNP (yellow stick model).

Notes: The residues are shown as wireframe colored by their element, and the hydrogen bonds are illustrated as dotted green lines. The water molecule in both (A) and (B) are shown as small red colored balls.

Figure 6 Docking of designed hybrid compounds at ATP binding site of P. aeruginosa GyrB. (A) Binding pose of the top scored compound PD_1a (grey stick model); (B) PD_1a and PD_1b (stick model colored by their element) superimposed with ADPNP (yellow stick model).

Notes: The residues are shown as wireframe colored by their element, and the hydrogen bonds are illustrated as dotted green lines. The water molecule in both (A) and (B) are shown as small red colored balls.

The 43 kDa E. coli GyrB subunit and the 43 kDa topoisomerase IV ParE subunit structures are strikingly similar in their conserved active sites, suggesting that inhibitors might simultaneously target both type II topoisomerase enzymes with a possibility for broad spectrum action.⁴⁸ Hence, we modeled the topoisomerase IV enzyme of P. aeruginosa using the E. coli template structure (1S16) (Figures S2 and S4), and we performed similar docking calculations of the designed hybrid compounds in the ADPNP bound active site cavity of ParE. The docked conformations obtained were similar to GyrB, indicating a highly conserved binding mode (Figures S6 and S7), and suggesting that the designed compounds also bind to the ParE subunit of Topoisomerase IV.

Docking of methotrexate and proposed hybrid compounds at the folate binding site of DHFR

The 2,4-diaminopteridine ring of methotrexate is situated into a deep hydrophobic pocket and interacts strongly with the P. aeruginosa DHFR through five major hydrogen bonds (Figure 4B). The N1 and the 2-amino group of methotrexate form a pair of hydrogen bonds with Asp30. The N1 is protonated and the carboxylate of the acidic residue Asp30 is ionized, facilitating a favorable ionic interaction.⁴⁹ The 4-amino group forms strong hydrogen bond interactions with the main chain carbonyl oxygen atom of Ile8, Ile104, and the hydroxyl group of Tyr110. The interaction pattern of diaminopteridine ring structures of all antifolates are highly preserved among the known DHFR complexes in various species.⁴⁹ In addition to the hydrogen bonding, the diaminopteridine ring forms Van der Waals interactions with Phe34, Ile8, and Ile104. The residues Phe34, Leu53, and Leu57 provide Van der Waals contacts to the phenyl group in the p-amino benzoate ring of methotrexate. The carbonyl group of p-aminobenzoate forms a hydrogen bond with Arg55, and the glutamate tail of methotrexate lies on the protein surface with its α-carboxylic group forming strong salt bridge interactions with the positively charged Arg60 residue (Figure 4B). The docking results show that all the proposed hybrid compounds have a good fit into the cavity of the DHFR active site and are positioned in a way similar to the inhibitor (methotrexate). The major residues participating in hydrogen bond interactions with the designed inhibitors in the DHFR active site are tabulated in Table 4. The docked conformation of the top ranked hybrid compounds, PD_2a and PD_1a, are shown in Figures 7A and 8A. The hybrid compounds PD_2a and 1a docks at the active site cavity of P. aeruginosa DHFR with its hydroxyl substituted aromatic rings being held in the interior of a deep cleft, while the sulfonamide moiety of sulfamethoxazole extends out towards the entrance of the binding pocket forming hydrogen bonding and hydrophobic interactions. The hydroxyl groups in the phenyl ring of phytochemicals are involved in hydrogen bonding interactions with six residues located interior of the active site Asp30, Thr123, Ile8, Ala9, Tyr110, and Ile104, while the ring contacts residues Phe34, Ile8, Ile104, and Leu31 by Van der Waals interactions. The phenyl ring of sulfamethoxazole adopts a similar orientation to that of the phenyl ring of methotrexate and forms Van der Waals contacts with Phe34 Leu31, Leu53, and Leu57 residues, which results in good ligand-enzyme affinity. The sulfonyl group of sulfamethoxazole forms hydrogen bond acceptor contacts with Arg55 by overlaying the carbonyl group of the p-aminobenzoate moiety in methotrexate (Figures 7A and 8A). The oxazole ring of the sulfamethoxazole demonstrates hydrogen bond acceptor interactions with the guanidium group of Arg60 at the entrance of the pocket in an orientation similar to that of the α-carboxylic group in methotrexate. In the case of PD_1b and 2b (Figures 7B and 8B), they are positioned in a similar way to that of PD_1a and 2a and also maintain similar interactions with the residues located in the interior of the active site. The phenyl ring of phytochemicals in both PD_1b and 2b are well adjusted in the hydrophobic pocket, situated deep inside the active site cavity, while the benzene sulfonamide ring adopts an orientation extended towards the Phe34 residue forming π-π contact and thereby disallowing the hydrogen bond interactions between the sulfonyl group and the Arg55 residue. On the other hand, the pyrimidine ring of sulfadiazine moieties in PD_1b and 2b establishes hydrogen bond interactions with both Arg55 and Arg60 residues.

Figure 7 Docking of designed hybrid compounds at folate binding site of P. aeruginosa DHFR. (A) Binding pose of the top scored compound PD_2a (grey stick model); (B) PD_2a and PD_2b (stick model colored by their element) superimposed with methotrexate (yellow stick model), and cofactor NADP is represented as aqua stick model.

Note: The residues are shown as wireframe colored by their element, and the hydrogen bonds are illustrated as dotted green lines.

Figure 8 Docking of designed hybrid compounds at folate binding site of P. aeruginosa DHFR. (A) Binding pose of the top scored compound PD_1a (grey stick model); (B) PD_1a and PD_1b (stick model colored by their element) superimposed with methotrexate (yellow stick model), and cofactor NADP is represented as aqua stick model.

Note: The residues are shown as wireframe colored by their element, and the hydrogen bonds are illustrated as dotted green lines.

Molecular dynamics simulations

To evaluate the stability of the unbound enzyme complexes (I–GyrB, IV–DHFR and VII–ParE), ligand-bound enzyme complexes (II–ADPNP-GyrB, V–MTX-NADP-DHFR and VIII–ADPNP-ParE) and top ranked hybrid compound-bound enzyme complexes (III–PD_2a-GyrB, VI–PD_2a-NADP-DHFR and IX–PD_2a-GyrB) were analyzed for the conformational changes with respect to the initial structure by RMSD during 3000 ps simulation time (Figure 9A and B and Figure S8). The average RMSD values obtained for complex I, IV and VII were approximately 2.060 Å, 1.97 Å, and 2.001 Å, respectively; for complex II, V and VIII they were 1.924 Å, 2.095 Å, and 2.181 Å, respectively; and for complex III, VI and XI they were 2.087 Å, 1.854 Å and 2.113 Å, respectively. As can be seen in the plots, after 500 ps, all the complexes tend to converge with low values of deviation of the molecules with respect to the original atomic positions, indicating the systems are stable and equilibrated. Furthermore, analyses of the RMSF of key active site residues versus the residue number for all the complexes are shown in Figure 10A–D. The RMSF value illustrates the average displacement of each residue in relation to their average backbone structure over the whole simulation and indicates the relative flexibility that describes each of these residue elements. As expected, the RMSF plots of all the complexes indicate that the main backbone residues with high fluctuations are those in the flexible loops. The residue positions between 78 and 88, and the turn which includes residues in positions between 248 and 251 of complexes I–GyrB, II– ADPNP-GyrB, and III–PD_2a-GyrB shows high fluctuations (Figure 10A). Similarly, the complexes IV–DHFR, V–MTX-NADP-DHFR, and VI–PD_2a-NADP-DHFR show loops that include residues between 67 and 78, and the α-helix structure which includes residues in positions between 90 and 96, which is located away from the binding site cavity, show high fluctuations (Figure 10C). On the other hand, the complexes VII–ParE, VIII–ADPNP-ParE, and IX–PD_2a-ParE indicate similar RMSF distributions and dynamic features, except in coils located between residue positions 71 and 81, 93 and 103, and turns 233–236, 299 and 303, which are further away from the active binding site areas (Figure S9). Both the systems with complexes I–III and VII–IX correspond to similar high fluctuations, and closer visual inspection reveals that conformational changes, due to the binding of respective similar substrates and inhibitors, produce displacement of the flanking loops/coils and turns, away from the binding site. The analysis of the mobile flexibilities of key active site residues RMSF values of active site residues in each complex fluctuates widely from 0.5 Å to 2.4 Å, indicating that the active site residue elements maintain their positions in relation to their corresponding enzyme complexes (Figures 10B, D, and S9). The comparison between the initial unbound and ligand-bound and inhibitor-bound complexes clearly shows the conformational changes favoring the binding of the substrates designed hybrid compounds in the active site cavities of the enzyme complexes. Thus, the low RMSD and RMSF deviations indicate that the inhibitor-enzyme complexes are stable and preserve the inhibitor interactions shown in the molecular docking experiments.

Figure 9 Root mean square deviation (RMSD) of backbone C^α atoms of the complexes versus time in picoseconds. (A) Unbound GyrB complex – I, ADPNP bound GyrB complex – II and top scored hybrid compound PD_2a bound GyrB complex-III; (B) Unbound DHFR complex – IV, methotrexate and cofactor NADP bound DHFR complex – V and top scored hybrid compound PD_2a and cofactor NADP bound DHFR complex – VI.

Figure 10 Root mean square fluctuation (RMSF) of backbone C^α atoms of the complexes versus residue number in the sequence. (A) Unbound GyrB complex – I, ADPNP bound GyrB complex – II and top-scored hybrid compound PD_2a bound GyrB complex – III; (B) Unbound DHFR complex – IV, methotrexate and cofactor NADP bound DHFR complex – V and top-scored hybrid compound PD_2a and cofactor NADP bound DHFR complex – VI; (C) Closer look of GyrB active site residues RMSF during the simulation of three complexes; (D) Closer look of DHFR active site residues RMSF during the simulation of three complexes.

Discussion

As a first step in the multitarget therapeutics design, the most crucial part is to identify or prioritize possible combination of targets that can be targeted by a single agent to achieve optimal therapeutic benefit.⁸ The network-based knowledge from scientific literature is increasingly gaining interest for the insights it offers into identification of relevant target sets and development of multitarget therapeutics.⁸ The knowledge-based methods employ text mining of biological data from the scientific literature ie, exploit information from literature to generate novel hypotheses, and connecting concepts by new and indirect associative techniques to contribute to drug discovery.^50,51 In light of the meager antibiotic pipeline for gram-negative bacteria and the threat of MDR, the motivation for rational design of multitarget inhibitors came from one of the best studied and investigated plant polyphenol, epigallocatechin gallate (from green tea). Epigallocatechin gallate and its derivatives have numerous biological activities including interactions with multiple enzymes, such as GyrB/ParE DHFR.^30–32 It is reported that ligands with exact and similar structures are more likely to have similar properties and therefore often bind proteins with similar active sites, and based on this, proteins can be connected for multitarget drug design.^52,53 Previously, using drug–phytochemical combination assays, we proposed that phytochemicals (protocatechuic acid, gallic acid, quercetin and myricetin), which are structural analogs and fragments of epigallocatechin gallate, bind to the folate binding site of DHFR and ATP-binding site of GyrB/ParE enzymes of P. aeruginosa and are synergistic in combination with sulfonamides, such as sulfamethoxazole and sulfadiazine.^33,34,54 More importantly, the major concern in our study was achieving higher minimum inhibitory concentration (MIC) values of the chosen phytochemicals compared to that of the antibiotics used. In support of this, several studies have reported that phytochemicals have less stability and selectivity, adverse phytochemical–antibiotic interactions, less bioavailability, and are less potent due to their MIC ranges (<1,000 μg/mL), in comparison to the compounds from microbial origin (antibiotics).^55,56 Interestingly, phytochemicals can be converted to pharmacologically acceptable antimicrobial agents by optimization through structural modifications.⁵⁷

Natural products have long been a source of highly diverse lead structures for the drug development process, representing almost half of the drugs currently in clinical use.⁵⁸ The existing natural antibiotics are either derived directly or have been successfully modified into broad spectrum compounds such as ampicillin, derived from penicillin and azithromycin derived from erythromycin.⁵⁹ Similarly, structural modification of plant-derived antibacterials may help to develop more potent, target-specific, pharmacologically acceptable, broad spectrum compounds.^57,59 In continuation of our previous herb–drug combination studies, the present study used 3D structure-based pharmacophore information for two different enzymes (GyrB/ParE and DHFR) from published literature, coupled with our previous phytochemical–drug combination assay data. We designed and validated a novel class of hybrid compounds (“phyto-drug”) for simultaneous inhibition of these enzymes using various in silico methods.

The target enzymes GyrB/ParE and DHFR are independently essential from two pharmacologically important pathways in the biological network whose perturbation could result in the desired therapeutic outcome. DNA gyrase and topoisomerase IV are type II topoisomerase enzymes involved in bacterial DNA replication and cell division, and their structural similarity typically allows ligands to inhibit both the enzymes.³² Despite being highly homologous enzymes, DNA gyrase composed of two heterodimeric subunits GyrA/GyrB primarily acts in initiation of DNA replication and elongation of nascent DNA, while topoisomerase IV (also consisting of two heterodimeric subunits ParC/ParE) functions to relax DNA supercoils and decatenates daughter chromosomal DNA.⁴⁷ Aminocoumarins, novobiocin and cyclothialidines, are the known inhibitors of ATPase catalytic region in the N-terminal domains of both GyrB/ParE. However, novobiocin has been withdrawn from clinical use due to toxicity reasons, and cyclothialidines are reported to have weak antibacterial activity due to poor penetration effects.⁴⁷ Similarly, DHFR is another important enzyme from the folic acid biosynthesis pathway that is involved in the reduction of dihydrofolate to tetrahydrofolate and is inhibited by diaminopyrimidine antibiotics such as trimethoprim and methotrexate.⁶⁰ Currently, the standard course of treatment for uncomplicated forms of urinary tract infections as well as prophylaxis and Pneumocystis carinii pneumonia is a combination of trimethoprim and sulfonamide drugs resulting in a synergistic action.⁶¹ Sulfonamides act as structural analogs of p-aminobenzoic acid, competitively inhibiting the additional enzyme dihydropteroate synthase in the folate pathway, which catalyzes the addition of p-aminobenzoic acid to dihydropterin pyrophosphate to form pteroic acid.⁶¹ Due to the low outer membrane permeability, poor affinity, and characteristic intrinsic resistance by efflux pumps, the antibacterial potency of trimethoprim and sulfonamides are limited (with MICs usually in the resistance range).⁶² Thus, the rational design and development of novel hybrid entities that can inhibit both GyrB/ParE and DHFR enzymes would be beneficial in treating various multidrug resistant phenotypes.

The designed hybrid compounds or “phyto-drug” evaluated for its physicochemical, drug likeness, and toxicity risk profiles using in silico-based calculations resulted in nontoxic behavior, highly drug-like with desirable physicochemical properties in comparison to that of individual phytochemicals and antibiotics (see Table 1). In addition, the structural factors involved in designed hybrids–receptor interactions evaluated using molecular orbitals and MEP profiles, based on stereoelectronic propertie appear to be in good agreement with the key pharmacophoric features required for multi-inhibition. An important aspect in our hybrid design framework is the approach for “revival of older antibiotics” and “molecular tinkering of phytochemicals” through structural combination design by linking common inhibitors (phytochemicals, protocatechuic acid and gallic acid) of DHFR and GyrB/ParE with synergistic sulfonamide antibiotics, which may help improve the properties and also bypass the mechanism of resistance of these compounds. Plant-derived phenolic acids (phytochemicals) have the potential to modify bacterial resistance by evading the efflux mechanisms in multidrug resistant gramnegative bacteria.⁶³ Also, the semisynthetic tailoring and repositioning of known classes of antibiotics will not only increase the efficacy by expanding the spectrum of use but will also offer incremental improvements over existing therapies in reducing resistance and toxicity to antibiotics.⁶⁴ Another, novel strategy used here is the integration of multispecificity in enzymes of unrelated pathways, as they are highly advantageous in overcoming resistance of targets in different pathways by significantly delaying resistance development. The binding modes of the designed hybrids at ATP-binding sites of GyrB/ParE and the folate-binding site of DHFR were calculated using docking and dynamics simulation and were shown to have conserved and stable binding affinity maintaining key molecular interactions necessary for multienzyme inhibition. These compounds are potential drug candidates, and hence identification and modification of the most biologically active promiscuous plant-derived natural products/antibiotics with the view to improve and balance their drug-like properties will offer novel agents with superior therapeutic properties against infectious diseases.

Conclusion

In summary, our findings, for the first time, have demonstrated the integration of knowledge-based methods, phytochemical–antibiotic interaction network data, and a fragment combination approach in the design and development of novel hybrid compounds, based on the conjugation of phytochemicals and antibiotics for blockage of enzymes (GyrB/ParE and DHFR), in order to restrain the MDR capability of P. aeruginosa. The biological rationale of our hybrid design framework is that it is much more difficult for a pathogen to develop resistance when an inhibitor has activity against multiple targets in two pharmacologically important pathways coupled with restriction of unwanted, off-target effects of the phytochemicals, due to its conjugation with antibiotics through structural modifications. It is also significant that our design strategy has opened an opportunity for revival of older antibiotics to bypass the resistance mechanisms through fragment combination. Finally, the identification of targets with common phytochemicals inhibitory profile through systemic interaction patterns, and network analysis through phytochemical–antibiotic combination assays may help to cluster common drug targets and form a base for accelerating the process of antibacterial multitarget drug discovery.

Supplementary materials

Sequence analysis and homology modeling

The amino acid sequences of the target enzymes P. aeruginosa DNA Gyrase subunit B (GyrB, Locus ID – PA0004), topoisomerase subunit B (ParE, Locus ID – PA4967) and dihydrofolate reductase (FolA, Locus ID – PA0350) were retrieved from Pseudomonas genome database (http://www.pseudomonas.com/). The amino acid sequences were subjected to BLASTP search against Protein Data Bank (PDB) and SWISS-MODEL template selection tool to identify the suitable templates for comparative modeling. Based on the template selection analysis, E. coli GyrB (PDB code: 1EI1, 2.3 Å), E. coli ParE (PDB code: 1S16, 2.1 Å) and E. coli DHFR (PDB code: 1RX3, 2.2 Å) were selected as suitable templates with sequence identity of 74%, 75%, and 45%, respectively. The alignment of P. aeruginosa GyrB/ParE and DHFR protein sequences and template structures – E. coli GyrB (1EI1), ParE (1S16), and DHFR (1RX3) protein sequences, were performed using CLUSTALW2 (http://www.ebi.ac.uk/Tools/clustalw2/index.html). The pair-wise alignment (Figure S1) was carefully examined for any potential alignment errors and the 3D model building of GyrB/ParE and DHFR was performed by two automated modeling programs: SWISS-MODEL¹ and MODELLER 9v7.² The generated homology models were submitted to WHATIF server (http://swift.cmbi.ru.nl/servers/html/index.html) to model missing side chains and remove atomic clashes (bumps).³ The RMSD values were obtained by structural superimposition of all the GyrB/ParE and DHFR homology models from both SWISS-MODEL and MODELLER with their respective E. coli templates (1EI1, 1S16, and 1RX3) using SUPERPOSE⁴ web server (http://wishart.biology.ualberta.ca/SuperPose/). Out of all the models estimated, the SWISS-MODEL constructed models resulted in low RMSD values (Figure S2), 0.46 Å (GyrB), 0.36 (ParE), and 0.39 Å (DHFR), respectively. RMSD values of ≤2 Å from the experimentally determined structures are considered as “highly accurate” models, while models with Cα RMSD above this threshold and ≤4 Å were termed “reliable” models, and such models are not fit for drug design and binding affinity analysis.⁵ Thus, low RMSD values obtained shows the selected homology models are acceptable and quite similar to their respective templates.

The quality of obtained P. aeruginosa GyrB/ParE and DHFR models were further subjected to evaluations based on stereochemical analysis, absolute measures on the geometry, and energy profiles independently using the structure assessment tools, PROCHECK, ANOLEA, and QMEAN available in the SWISS-MODEL Workspace,¹ Verify3D⁶ and ProSA.⁷ The Ramachandran plot produced from PROCHECK analysis (Figure S3, S4, and S5) showed 84.3% (GyrB), 90.5% (ParE), and 86.5% (DHFR) amino acid residues within the most favored regions, 15.7% (GyrB), 9.5% (ParE), and 12.7% (DHFR) within additionally allowed regions indicating that the models have good overall stereochemical quality. The QMEAN analysis used for the estimation of the global quality of the entire model as well as for the local per residue analysis of different regions within the models resulted in the QMEAN6 global score values 0.707 (GyrB), 0.598 (ParE), and 0.752 (DHFR), which are well within the quality range (0–1) with higher values for better models.⁸ In addition, estimated QMEAN Z-scores represent a measure of the absolute quality of the models by providing an estimate of scores obtained for high-resolution reference structures solved experimentally by X-ray crystallography models. Low quality are expected to have strong negative Z-scores, and “good structures” have (mean Z-scores = −0.65).⁹ The observed Z-score values for the predicted models were −0.718 (GyrB), −0.53 (ParE), and −0.077 (DHFR), suggesting these models are of comparable quality to high-resolution experimental structures. The energy profiles of the models were estimated using both ANOLEA and ProSA analysis showing regions with predominant low energy (negative values). The ProSA tool is widely used to check 3D models of protein structures for potential errors by calculating z-score and a plot of its residue energies, z-score value, is displayed in a plot that contains the z-scores of all experimentally determined protein chains in current PDB and models with z-scores outside the range characteristic for native proteins indicate erroneous structures, while the energy plot shows the local model quality by plotting residue energies as a function of amino acid sequence position in which positive values correspond to problematic or erroneous parts of a model.⁷ The calculated z-scores −8.45 (GyrB), −8.7 (ParE), and −5.19 (DHFR) of the modeled structures are within the range of scores typically found for native proteins of similar size and the energy plot showed all the residues in GyrB/ParE have negative energies, while residues in DHFR showed maximum residues with negative energy with very few residues displaying positive energy values (Figure S3, S4, and S5). Finally, Verify3D was used to test the accuracy of a 3D model by comparing the model to its own amino acid sequence using a 3D profile.⁶ The average 3D-1D profile scores for the model (GyrB, 0.67) and that of its template (1EI1, 0.71), model (ParE, 0.73) and its template (1S16, 0.71) similarly, (DHFR, 0.69) and its template (1RX3, 0.71) are very close. Hence, the validation based on the geometry, energy profiles and overall stereo-chemical quality suggests that the homology models obtained are appropriate for further protein-ligand interaction studies.

Figure S1 Alignment of target (P. aeruginosa) sequence and E. coli template protein sequence. (A) DNA Gyrase subunit B (GyrB) with template (PDB: 1EI1). (B) Topoisomerase IV subunit B (ParE) with template (PDB: 1S16). (C) Dihydrofolate reductase (DHFR) with template (PDB: 1RX3).

Figure S2 Superimposition between P. aeruginosa homology models and the template crystal structures. (A) GyrB model (blue) and template structure E. coli (1EI1, red) bound with ADPNP (yellow stick model). (B) ParE model (blue) and template structure E. coli (1S16, green) bound with ADPNP (magenta stick model). (C) DHFR model (blue) and template structure E. coli (1RX3, green) bound with the inhibitor methotrexate (yellow stick model), and NADP cofactor (aqua stick model).

Note: Active site residues are shown as wireframe colored by their element.

Figure S3 Validation of P. aeruginosa GyrB homology model (A) Ramachandran plot and (B) ProSA (protein structure analysis) energy plot.

Figure S4 Validation of P. aeruginosa ParE homology model (A) Ramachandran plot and (B) ProSA (protein structure analysis) energy plot.

Figure S5 Validation of P. aeruginosa DHFR homology model (A) Ramachandran plot and (B) ProSA (protein structure analysis) energy plot.

Figure S6 Docking of designed hybrid compounds at ATP binding site of P. aeruginosa ParE. (A) Binding pose of the top-scored compound PD_2a (grey stick model) (B) PD_2a and PD_2b (stick model colored by their element) superimposed with ADPNP (magenta stick model).

Notes: The residues are shown as wireframe colored by their element and labeled in white. The hydrogen bonds are illustrated as dotted green lines. The water molecule in both (A) and (B) is shown as a red colored ball model.

Figure S7 Docking of designed hybrid compounds at ATP binding site of P. aeruginosa ParE. (A) Binding pose of the top-scored compound PD_1a (grey stick model). (B) PD_1a and PD_1b (stick model colored by their element) superimposed with ADPNP (magenta stick model).

Notes: The residues are shown as wireframe colored by their element and labeled in white. The hydrogen bonds are illustrated as dotted green lines. The water molecule in both (A) and (B) is shown as a red colored ball model.

Figure S8 Root mean square deviation (RMSD) of backbone C^α atoms of the complexes versus time in picoseconds. Unbound ParE complex – VII, ADPNP bound ParE complex – VIII and top-scored hybrid compound PD_2a bound ParE complex – IX.

Figure S9 Root mean square fluctuation (RMSF) of backbone C^α atoms of the complexes versus residue number in the sequence. (A) Unbound ParE complex – VII, ADPNP bound ParE complex – VIII and top-scored hybrid compound PD_2a bound ParE complex – IX; (B) Closer look of ParE active site residues RMSF during the simulation of three complexes.

Table S1 Docking results of ADPNP and designed hybrid compounds (ranked on the basis of binding affinity scores) at the ATP binding site in P. aeruginosa ParE

Compounds	Scoring function		EInter^c	EIntra^d	Hbond^e	LE1^f	LE3^g	Binding affinity (kJ/mol)	Residues involved in hydrogen bonding
	MolDock^a	Rerank^b
PD_2a	−190.47	−133.7	−200.73	10.25	−25.74	−6.35	−4.46	−37.18	Asp68, Ser42, Glu45, Tyr104, Asn41, Ile115, Gly114, Gly112, His111, Leu110, HOH
PD_1a	−193.98	−152.5	−204.23	10.25	−21.81	−6.26	−4.92	−33.83	Asp68, Ser42, Glu45, Tyr104, Asn41, Ile115, Gly114, Gly112, His111, Leu110, HOH
PD_2b	−171.55	−129.94	−191.33	19.78	−19.58	−5.72	−4.33	−33.01	Asp68, Ser42, Glu45, Tyr104, Asn41, Ile115, Gly114, Gly112, His111, Leu110, HOH
PD_1b	−178.9	−150.87	−197.23	18.33	−18.37	−5.77	−4.87	−31.3	Asp68, Ser42, Glu45, Tyr104, Asn41, Ile115, Gly114, Gly112, His111, Leu110, Val113, HOH
ADPNP	−234.21	−184.57	−224.59	−9.62	−32.15	−7.56	−5.95	−33.86	Asp68, Tyr104, Gly97, Lys98, Asn41, Ile115, Gly114, Val113, Gly112, His111, Leu110, Lys333, HOH

Abbreviations: MolDock, docking score function by MolDock; Rerank, reranking score function; E_Inter, ligand – protein (GyrB) interaction energy; E_Intra, internal energy of the ligand; HBond, hydrogen bonding energy; LE1, Ligand efficiency 1 – MolDock score divided by Heavy atoms count; LE3, Ligand efficiency 3 – Rerank score divided by Heavy atoms count; ADPNP, (nonhydrolyzable analog of ATP) is GyrB substrate.

References

• Arnold K Bordoli L Kopp J Schwede T The SWISS-MODEL workspace: a web-based environment for protein structure homology modelling Bioinformatics 2006 22 2 195 201 16301204
   PubMed Web of Science ®Google Scholar
• Sali A Blundell TL Comparative protein modelling by satisfaction of spatial restraints J Mol Biol 1993 234 3 779 815 8254673
   PubMed Web of Science ®Google Scholar
• Rodriguez R Chinea G Lopez N Pons T Vriend G Homology modeling, model and software evaluation: three related resources Bioinformatics 1998 14 6 523 528 9694991
   PubMed Web of Science ®Google Scholar
• Maiti R Van Domselaar GH Zhang H Wishart DS SuperPose: a simple server for sophisticated structural superposition Nucleic Acids Res 2004 32 Web Server issue W590 W594 15215457
   PubMed Web of Science ®Google Scholar
• Rayan A New tips for structure prediction by comparative modeling Bioinformation 2009 3 6 263 267 19255646
   PubMedGoogle Scholar
• Eisenberg D Lüthy R Bowie JU VERIFY3D: assessment of protein models with three-dimensional profiles Methods Enzymol 1997 277 396 406 9379925
   PubMed Web of Science ®Google Scholar
• Wiederstein M Sippl MJ ProSA-web: interactive web service for the recognition of errors in three-dimensional structures of proteins Nucleic Acids Res 2007 35 Web Server issue W407 W410 17517781
   PubMed Web of Science ®Google Scholar
• Benkert P Tosatto SC Schomburg D QMEAN: a comprehensive scoring function for model quality assessment Proteins 2008 71 1 261 277 17932912
   PubMed Web of Science ®Google Scholar
• Benkert P Biasini M Schwede T Toward the estimation of the absolute quality of individual protein structure models Bioinformatics 2011 27 3 343 350 21134891
   PubMed Web of Science ®Google Scholar

Acknowledgments

Premkumar Jayaraman, Kishore R Sakharkar, and Meena K Sakharkar conceived the study, designed the methods, carried out the data analysis, and drafted the manuscript. Meena K Sakharkar supervised the research and contributed to drafting the manuscript, Lim ChuSing Daniel, Kishore R Sakharkar and Sarinder K Dhillon, and Mohammad Imran Siddiqi gave useful suggestions and improved on the manuscript. All authors have read and approved the final manuscript.

This work was funded by a seed grant from Biopharmaceutical Engineering Cluster, Nanyang Technological University, Singapore and an AcRF Tier1 grant from Nanyang Technological University, Singapore.

Disclosure

The authors report no conflicts of interest.