In Silico Design, Molecular Docking And ADME-Toxicity Evaluation of Substituted Benzimidazole Derivatives As Potential Antifungal Agents

Fungal infections have emerged as a critical area of unmet medical need, particularly among immunocompromised populations including organ transplant recipients, cancer patients undergoing chemotherapy, individuals with HIV/AIDS, and those on prolonged corticosteroid therapy [1]. The spectrum of clinically significant fungal pathogens encompasses species of Candida, Aspergillus, Cryptococcus, and dermatophytic fungi, collectively responsible for an estimated 1.5 million deaths annually worldwide [2]. Despite this considerable disease burden, the antifungal drug armamentarium remains limited to a handful of mechanistic classes, namely the azoles targeting ergosterol biosynthesis, polyenes disrupting membrane integrity, echinocandins inhibiting beta-glucan synthesis, and antimetabolites such as flucytosine [3]. This narrow therapeutic landscape is further constrained by intrinsic and acquired resistance mechanisms that have progressively eroded the clinical utility of frontline agents. Azole resistance in Candida glabrata and Candida auris, along with the emergence of multi-azole- resistant Aspergillus fumigatus strains, underscores the urgency for new antifungal scaffolds acting through alternative molecular mechanisms [4].

Tubulin, the structural protein component of microtubules, has long been recognised as a selective antifungal target, primarily because of meaningful differences in the colchicine-binding site between fungal and mammalian tubulin isoforms [5]. Microtubules are dynamic polymers formed by the polymerisation of alpha-beta tubulin heterodimers and are indispensable for mitotic spindle assembly, chromosome segregation, and intracellular transport [6]. Agents that bind the colchicine site stabilise the curved conformation of the tubulin heterodimer, thereby preventing longitudinal contacts required for polymerisation [7]. Benzimidazole antifungals such as mebendazole, albendazole, and thiabendazole exploit this mechanism, binding selectively to fungal tubulin owing to structural divergences at residues that form the colchicine pocket, thereby achieving antiparasitic and antifungal activity at doses that spare mammalian microtubules to a meaningful extent [8]. However, the pharmacokinetic limitations of classical benzimidazoles, including poor aqueous solubility, erratic gastrointestinal absorption, and systemic adverse effects, highlight the need for structurally optimised analogues.

Benzimidazole is a privileged pharmacophore in medicinal chemistry, incorporating a benzene ring fused to an imidazole nucleus to create a bicyclic aromatic system capable of extensive hydrogen bonding, aromatic interactions, and electronic delocalisation [9]. The scaffold has been exploited across a broad range of biological activities including antifungal, antiparasitic, antiviral, anticancer, and anti-inflammatory applications [10]. Structural modifications at the N-1 position, C-2 position, and through Schiff base condensation with the hydrazinyl function have been demonstrated to substantially modulate binding affinity, selectivity, and physicochemical properties [11]. The introduction of hydrazone linkages connects the benzimidazole core to a substituted terminal aromatic ring, creating a conjugated framework that increases planarity and enhances pi-stacking capacity within binding pockets [12].

Computational drug discovery approaches have transformed the efficiency and directionality of the lead identification phase [13]. Molecular docking provides a mechanistic understanding of ligand-protein complementarity at the atomic level, while in silico ADME prediction and toxicity profiling through tools such as SwissADME and ProTox-II allow for the rapid elimination of pharmacokinetically inferior or potentially toxic candidates at early stages of the discovery pipeline [14,15]. These computational filters are invaluable for prioritising compounds for synthesis and experimental validation, significantly reducing both cost and time-to-candidate [16].

Given these considerations, the present study was designed to rationally construct a series of eight benzimidazole-based Schiff base hydrazone derivatives carrying diverse substituents, including bromo, fluoro, iodo, nitro, trifluoromethyl, hydrazinyl, cyano, and carboxylic acid groups, on both the central phenyl aryl ring and the terminal aromatic moiety. These compounds were then evaluated computationally using the 1SA0 crystal structure of the tubulin-colchicine stathmin-like domain complex as the molecular target. The objectives were to determine binding affinities relative to mebendazole, analyse structure-activity relationships at the electronic and steric level, assess oral drug-likeness according to established pharmacokinetic rules, and evaluate the preliminary toxicological risk profile, collectively providing a rational computational basis for the advancement of optimised benzimidazole antifungal leads.

2. MATERIALS AND METHODS

2.1 Target Protein Selection and Preparation

The crystal structure of the tubulin-colchicine stathmin-like domain complex (PDB ID: 1SA0) was retrieved from the RCSB Protein Data Bank (https://www.rcsb.org). This structure, resolved at 3.58 Å by X-ray crystallography, contains the alpha-beta tubulin heterodimer bound to colchicine within the interdimer interface, representing the classical colchicine-binding site that is targeted by benzimidazole antifungals [7]. The protein was imported into BIOVIA Discovery Studio Visualizer 2021 for structural inspection, removal of crystallographic water molecules, co-crystallised ligands, and ionic cofactors not essential to binding pocket architecture. Polar hydrogen atoms were added, and Gasteiger partial charges were assigned to all protein atoms to ensure accurate electrostatic representation during the docking calculation. Energy minimisation of the prepared protein was performed to relieve steric clashes arising from crystal packing forces prior to grid box generation.

Figure 1. Three-dimensional structure of the tubulin-colchicine stathmin-like domain complex (PDB ID: 1SA0) used for molecular docking studies.

2.2 Ligand Design and Structure Preparation

The eight target molecules were designed as substituted benzimidazole-Schiff base hydrazone conjugates featuring a 2-arylbenzimidazole core coupled through a hydrazone linkage (C=N-NH) to a variably substituted terminal aryl group. Structural diversity was introduced by positioning different substituents on both the central aryl ring attached to the C-2 position of benzimidazole and on the terminal phenyl ring of the hydrazone moiety. Specifically, the central ring carried either a 4-methyl, 5-nitro, 5-trifluoromethyl, 5-hydrazinyl, or 5-carboxylic acid group, while the terminal ring featured 4-bromo, 3-bromo, 4-iodo, 3-fluoro, 4-fluoro, 3-chloro, 4-cyano, or 3-iodo substituents. Two-dimensional chemical structures were drawn using ACD/ChemSketch (Advanced Chemistry Development, Inc.), and canonical SMILES strings were generated for each compound. Structures were subsequently imported into PyRx 0.8 (Virtual Screening Software) where they were converted to three-dimensional format using the integrated Open Babel toolbox. Each ligand underwent geometry optimisation using the Universal Force Field (UFF), followed by energy minimisation to achieve the lowest-energy three-dimensional conformation, which was then converted to the PDBQT format required for AutoDock Vina docking.

Mebendazole (C₁₆H₁₃N₃O₃; MW 295.29 g/mol), the WHO-listed benzimidazole reference antifungal drug, was retrieved from the PubChem database and subjected to identical preparation procedures to serve as the positive control.

2.3 Molecular Docking Procedure

All docking simulations were carried out using AutoDock Vina as implemented in PyRx 0.8, which employs an iterated local search global optimiser combined with a knowledge-based empirical scoring function to estimate binding free energies in kcal/mol [13]. A grid box was defined around the colchicine-binding site of 1SA0 at coordinates corresponding to the centroid of the co-crystallised colchicine molecule. The grid dimensions were set at 25 x 25 x 25 Å with a spacing of 0.375 Å to adequately encompass the binding pocket and accommodate the planned molecular scaffolds. The exhaustiveness parameter was set to 8, generating nine binding modes per ligand from which the conformation displaying the lowest binding free energy was selected as the representative docked pose. All docking results were validated by re-docking colchicine into the binding site and verifying that the re-docked pose reproduced the crystallographic binding conformation within an acceptable root-mean-square deviation threshold.

2.4 Interaction Analysis

Detailed analysis of non-covalent interactions between each docked ligand and the binding site residues was performed using BIOVIA Discovery Studio Visualizer. Two-dimensional ligand-protein interaction diagrams were generated to identify hydrogen bond donors and acceptors, hydrophobic contacts, pi-pi stacking interactions, halogen bonds, and electrostatic interactions. Three-dimensional binding poses were visualised to assess the spatial orientation of each ligand within the active site cavity and to understand the structural basis for observed differences in binding affinity across the series.

2.5 ADME and Drug-Likeness Prediction

Pharmacokinetic and physicochemical properties were predicted for all eight designed molecules and mebendazole using the SwissADME web server (http://www.swissadme.ch) [14]. Input was provided as canonical SMILES strings. Parameters evaluated included molecular weight, topological polar surface area (TPSA), number of hydrogen bond acceptors and donors, count of rotatable bonds, lipophilicity (consensus log P), gastrointestinal (GI) absorption, blood-brain barrier (BBB) permeation, and ESOL-based aqueous solubility class. Drug-likeness was evaluated against the Lipinski Rule of Five [17], Veber criteria [18], and Egan filter [19]. Bioavailability score was calculated from the radar plot parameters provided by SwissADME. Medicinal chemistry filters including PAINS (Pan Assay Interference Compounds) and Brenk structural alerts were applied to identify potentially promiscuous or reactive substructures [14].

2.6 Toxicity Prediction

In silico toxicological profiling was conducted using the ProTox-II webserver (https://tox-new.charite.de/protox_II/) [15], which employs machine learning models trained on data from established toxicological databases to predict acute oral toxicity LDâ‚…â‚€ values and GHS toxicity classification. Additional endpoint predictions were obtained for hepatotoxicity, nephrotoxicity, neurotoxicity, respiratory toxicity, cardiotoxicity, carcinogenicity, immunotoxicity, mutagenicity, and cytotoxicity. Probability scores for each endpoint were recorded alongside qualitative active-inactive classifications to provide a comprehensive preliminary safety assessment for each compound.

3. RESULTS AND DISCUSSION

3.1 Molecular Design and Structural Characteristics

Figure 2. Chemical structures of the designed substituted benzimidazole derivatives (Molecules 1–8).

Figure 3. Two-dimensional interaction diagram of Molecule 1 within the active site of 1SA0.

Figure 4. Two-dimensional interaction diagram of Molecule 2 within the active site of 1SA0.

Figure 5. Two-dimensional interaction diagram of Molecule 3 showing key hydrogen bonding, halogen bonding, and π–π interactions within the colchicine-binding site of tubulin.

Figure 6. Two-dimensional interaction diagram of Molecule 4.

Figure 7. Two-dimensional interaction diagram of Molecule 5.

Figure 8.  Two-dimensional interaction diagram of Molecule 6.

Figure 9. Two-dimensional interaction diagram of Molecule 7.

Figure 10. Two-dimensional interaction diagram of Molecule 8 within the colchicine-binding site of 1SA0.

Figure 11. Two-dimensional interaction diagram of the reference drug mebendazole docked within the active site of 1SA0.

Figure 12.  Three-dimensional binding pose of Molecule 3 within the colchicine-binding pocket of tubulin (1SA0).

Figure 13. Three-dimensional binding pose of mebendazole within the active site of tubulin (1SA0).

Acute oral toxicity predictions from ProTox-II reveal a notable finding: seven of the eight designed molecules are assigned to GHS Toxicity Class 4, corresponding to predicted LDâ‚…â‚€ values in the range of 300 to 2000 mg/kg, while mebendazole is classified as Toxicity Class 2 with a predicted LDâ‚…â‚€ of only 35 mg/kg. This suggests that, computationally, the designed benzimidazole derivatives carry a considerably reduced acute oral toxicity risk compared with the reference standard, which is a meaningful differentiation in the context of the therapeutic window for antifungal agents. The exception is Molecule 4, predicted with a Toxicity Class 1 classification and an alarmingly low LDâ‚…â‚€ of 4 mg/kg, which warrants careful experimental evaluation to determine whether this extreme toxicity prediction reflects genuine structural liabilities associated with the trifluoromethyl-fluoro substitution pattern or represents a limitation of the computational model due to low structural similarity with the training set compounds (average similarity 34.73%, the lowest in the series).

Hepatotoxicity and neurotoxicity are predicted as active for all eight designed molecules, a finding that is consistent with the broad bioactivity profile of nitrogen-rich aromatic heterocycles and should be considered a high-priority endpoint for experimental safety assessment in any future in vitro or in vivo studies. In this respect, it is important to acknowledge that in silico toxicity prediction tools such as ProTox-II operate on machine learning classifiers trained predominantly on acute rodent data and chemical property-based models, and their predictions carry significant uncertainty, particularly when test compound structural similarity to training set molecules is limited. The probability scores for hepatotoxicity range from 0.57 to 0.67 across the series, with mebendazole showing a score of 0.59, indicating that none of the designed molecules substantially deviate from the reference standard in terms of this particular endpoint prediction.

The key differentiator in the toxicity profile is observed in carcinogenicity and mutagenicity assessment. Molecules 3, 1, 5, 8, and 4 are predicted as carcinogenicity-inactive, while Molecules 2, 6, and 7 show carcinogenicity-active classification. Molecule 2, which also shows active mutagenicity, bears a 5-nitro group on the central phenyl ring, a widely recognised genotoxic structural alert capable of generating reactive nitroso and hydroxylamine metabolites under enzymatic reduction [22]. This strongly argues against further development of Molecule 2 and reinforces the general medicinal chemistry principle that aromatic nitro groups should be avoided in drug candidates wherever possible. Molecule 7, classified as both carcinogenicity-active and mutagenicity-active, carries a 4-cyano terminal substituent; the mutagenicity prediction for this compound may relate to the potential metabolic hydrolysis of the nitrile to an amino group, generating a primary aromatic amine metabolite that can undergo oxidative activation to reactive electrophilic species.

Molecule 3 stands out in the toxicological evaluation as well, being predicted as carcinogenicity-inactive, immunotoxicity-inactive, mutagenicity-inactive, and cytotoxicity-inactive, while showing the lowest organ toxicity burden among the halogenated molecules. Combined with its superior binding affinity of -9.0 kcal/mol, high GI absorption prediction, full Lipinski compliance, and absence of nephrotoxicity and cardiotoxicity signals, Molecule 3 emerges as the most favourable overall candidate in this series, balancing potency predictions with an acceptable computational safety profile.

Nephrotoxicity is predicted as inactive for six of the eight designed molecules. Only Molecule 8 and mebendazole show active nephrotoxicity predictions. The 5-carboxylic acid group in Molecule 8 may contribute to renal accumulation risk given the well-established renal elimination route for organic acid compounds. Cardiotoxicity is predicted as inactive across all designed molecules and for mebendazole, which is a reassuring finding given the significant concern around hERG potassium channel inhibition and cardiac arrhythmia risk that surrounds many nitrogen-containing aromatic scaffolds.

3.7 Comparative Assessment of Lead Molecule 3 Against Mebendazole

Viewed collectively, Molecule 3 demonstrates computational superiority over the reference standard mebendazole across several key parameters. Its binding affinity of -9.0 kcal/mol exceeds that of mebendazole by 0.3 kcal/mol, and the interaction profile shows a complementary combination of hydrogen bonding, hydrophobic contacts, pi-pi stacking with Phe255, and a unique halogen bond mediated by the 4-iodo group. High GI absorption is predicted for Molecule 3 compared to the low absorption of mebendazole, representing a potential pharmacokinetic advantage. Both molecules comply with the Lipinski Rule of Five. While mebendazole shows no PAINS or Brenk alerts, Molecule 3 triggers two PAINS and two Brenk flags, primarily reflecting the hydrazone bridge and the iodo substituent. From a toxicity perspective, Molecule 3 is assigned to the less toxic Class 4 compared to mebendazole's Class 2, and it avoids the carcinogenicity and mutagenicity flags that characterise the reference standard. These collective computational attributes establish Molecule 3 as a compelling lead for synthetic elaboration and experimental pharmacological evaluation.

4. FUTURE PERSPECTIVES

The findings of this computational study open several promising directions for continued investigation. The first priority is the chemical synthesis of the highest-ranked compounds, particularly Molecule 3, using established condensation protocols for benzimidazole-hydrazone formation, followed by spectroscopic characterisation. Experimental validation through minimum inhibitory concentration assays against a panel of clinically relevant fungal pathogens including Candida albicans, Candida glabrata, Aspergillus fumigatus, and dermatophytic species is essential to determine actual antifungal potency. Concurrent in vitro tubulin polymerisation inhibition assays would confirm the predicted mechanistic basis of activity. Cytotoxicity evaluation against mammalian cell lines, particularly hepatic and neuronal cells in light of the computational hepatotoxicity and neurotoxicity flags, would establish the therapeutic index and guide further structural optimisation.

From a structural optimisation perspective, the clear halogen bond advantage conferred by the iodo group may be further exploited through the introduction of non-classical halogen bond donors or by bioisosteric replacement of problematic fragments identified through PAINS alerts. Prodrug strategies or pharmaceutical formulation approaches, including nanoparticulate or cyclodextrin complexation, may be explored to address the poor aqueous solubility that affects the majority of compounds in this series. Molecular dynamics simulations over extended timescales would provide valuable insights into binding pose stability and the conformational dynamics of the protein-ligand complex, complementing the static snapshot provided by docking. Finally, expansion of the compound library by incorporating additional substitution patterns, particularly those combining the favourable 4-iodo terminal group with electron-donating groups at the 5-position of the central aryl ring, may yield analogues with further improved affinity and pharmacokinetic profiles.

CONCLUSION

The present computational investigation describes the rational in silico design and comprehensive evaluation of eight novel benzimidazole-Schiff base hydrazone derivatives as potential tubulin-targeting antifungal agents against the colchicine-binding site of the 1SA0 crystal structure. The systematic variation of substituents across the central and terminal aryl rings revealed clear and interpretable structure-activity relationships. Among the designed series, Molecule 3, bearing a 4-iodo terminal substituent and a 4-methyl central ring group, achieved the highest binding affinity of -9.0 kcal/mol, surpassing the reference antifungal mebendazole (-8.7 kcal/mol). The superior performance of Molecule 3 was attributed to a multivalent interaction profile incorporating hydrogen bonding with Lys352 and Val236, hydrophobic contacts with Ala314, pi-pi stacking with Phe255, and a directional halogen bond contributed by the para-iodo group. Pharmacokinetic profiling confirmed Lipinski compliance, high predicted GI absorption, and an acceptable bioavailability score of 0.55 for Molecule 3, while the ProTox-II toxicity analysis assigned it to the less hazardous GHS Class 4, with no carcinogenicity, mutagenicity, nephrotoxicity, cardiotoxicity, or cytotoxicity signals. The nitro-containing Molecule 2 was identified as the least suitable candidate due to its genotoxic alert profile, while Molecule 4 raised acute toxicity concerns requiring experimental clarification. Collectively, these computational findings provide a strong and scientifically coherent basis for the prioritisation of Molecule 3 as a lead compound for synthetic preparation and experimental antifungal evaluation, with in vitro minimum inhibitory concentration determination, tubulin polymerisation inhibition assays, and cell-based cytotoxicity profiling constituting the logical next experimental steps.

REFERENCES

1. Bongomin F, Gago S, Oladele RO, Denning DW. Global and multi-national prevalence of fungal diseases -- estimate precision. J Fungi (Basel). 2017;3(4):57. https://doi.org/10.3390/jof3040057
2. Bhattacharya S, Sae-Tia S, Fries BC. Candidiasis and mechanisms of antifungal resistance. Antibiotics. 2020;9(6):312. https://doi.org/10.3390/antibiotics9060312
3. Perfect JR. The antifungal pipeline: a reality check. Nat Rev Drug Discov. 2017;16(9):603-616. https://doi.org/10.1038/nrd.2017.46
4. Fisher MC, Hawkins NJ, Sanglard D, Gurr SJ. Worldwide emergence of resistance to antifungal drugs challenges human health and food security. Science. 2018;360(6390):739-742. https://doi.org/10.1126/science.aap7999
5. Blume-Jensen P, Hunter T. Oncogenic kinase signalling. Nature. 2001;411(6835):355-365. https://doi.org/10.1038/35077225
6. Jordan MA, Wilson L. Microtubules as a target for anticancer drugs. Nat Rev Cancer. 2004;4(4):253-265. https://doi.org/10.1038/nrc1317
7. Ravelli RBG, Gigant B, Curmi PA, et al. Insight into tubulin regulation from a complex with colchicine and a stathmin-like domain. Nature. 2004;428(6979):198-202. https://doi.org/10.1038/nature02393
8. Lacey E. The role of the cytoskeletal protein, tubulin, in the mode of action and mechanism of drug resistance to benzimidazoles. Int J Parasitol. 1988;18(7):885-936. https://doi.org/10.1016/0020-7519(88)90118-2
9. Keri RS, Hiremathad A, Budagumpi S, Nagaraja BM. Comprehensive review in current developments of benzimidazole-based medicinal chemistry. Chem Biol Drug Des. 2015;86(1):19-65. https://doi.org/10.1111/cbdd.12462
10. Ansari A, Ali A, Asif M, Shamsuzzaman. Review: biologically active pyrazole derivatives. New J Chem. 2017;41(1):16-41. https://doi.org/10.1039/C6NJ03245C
11. Sondhi SM, Rani R, Roy P, Agrawal SK, Saxena AK. Microwave-assisted synthesis, anticancer, and anti-inflammatory activity evaluation of Schiff's base and acetohydrazide derivatives of indole. Mol Divers. 2009;13(3):357-366. https://doi.org/10.1007/s11030-009-9126-8
12. Narasimhan B, Sharma D, Kumar P. Biological importance of the indole nucleus in recent years: a comprehensive review. Med Chem Res. 2012;21(1):1-21. https://doi.org/10.1007/s00044-010-9578-9
13. Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem. 2010;31(2):455-461. https://doi.org/10.1002/jcc.21334
14. Daina A, Michielin O, Zoete V. SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci Rep. 2017;7:42717. https://doi.org/10.1038/srep42717
15. Banerjee P, Eckert AO, Schrey AK, Preissner R. ProTox-II: a webserver for the prediction of toxicity of chemicals. Nucleic Acids Res. 2018;46(W1):W257-W263. https://doi.org/10.1093/nar/gky318
16. Leelananda SP, Lindert S. Computational methods in drug discovery. Beilstein J Org Chem. 2016;12:2694-2718. https://doi.org/10.3762/bjoc.12.267
17. Lipinski CA, Lombardo F, Dominy BW, Feeney PJ. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv Drug Deliv Rev. 2001;46(1-3):3-26. https://doi.org/10.1016/S0169-409X(00)00129-0
18. Veber DF, Johnson SR, Cheng HY, et al. Molecular properties that influence the oral bioavailability of drug candidates. J Med Chem. 2002;45(12):2615-2623. https://doi.org/10.1021/jm020017n
19. Egan WJ, Merz KM Jr, Baldwin JJ. Prediction of drug absorption using multivariate statistics. J Med Chem. 2000;43(21):3867-3877. https://doi.org/10.1021/jm000292e
20. Cavallo G, Metrangolo P, Milani R, et al. The halogen bond. Chem Rev. 2016;116(4):2478-2601. https://doi.org/10.1021/acs.chemrev.5b00484
21. Müller K, Faeh C, Diederich F. Fluorine in pharmaceuticals: looking beyond intuition. Science. 2007;317(5846):1881-1886. https://doi.org/10.1126/science.1131943
22. Purohit V, Basu AK. Mutagenicity of nitroaromatic compounds. Chem Res Toxicol. 2000;13(8):673-692. https://doi.org/10.1021/tx000002x