Theoretical activity prediction, structure-based design, molecular docking and pharmacokinetic studies of some maleimides against Leishmania donovani for the treatment of leishmaniasis

Background: Leishmaniasis is a neglected tropical disease caused by a group of protozoan of the genus Leishmani a and transmitted to humans majorly through the bite of the female sand fly. It is prevalent in the tropical regions of the world especially in Africa and estimated to affect a population of about 12 million people annually. This theoretical study was therefore conducted in support of the search for more effective drug candidates for the treatment of leishmaniasis. This study focuses on the in silico activity prediction of twenty-eight (28) maleimides, structure-based design, molecular docking study and pharmacokinetics analysis of the newly designed maleimides. All the studied compounds were drawn using ChemDraw Ultra and optimized by the density functional theory (DFT) approach using B3LYP with 6-31G⁄ basis set. Results: The built QSAR model was found to satisfy the requirement of both internal and external validation tests for an acceptable QSAR model with R 2 = 0.801, R 2adj = 0.748, Q 2cv = 0.710, R 2test = 0.892 and c R p2 = 0.664 and has shown excellent prediction of the studied compounds. Among the five (5) protein receptors utilized for the virtual docking screening, pyridoxal kinase (PdxK) receptor (Pdb id = 6k91) showed the strongest binding interactions with compounds 14 , 21 and 24 with the highest binding affinities of − 7.7, − 7.7 and − 7.8 kcal/mol, respectively. The selected templates ( 14 , 21 and 24 ) were used to design twelve (12) new compounds (N1–N12) with higher docking scores than the templates. N7 (affinity = − 8.9 kcal/mol) and N12 ( − 8.5 kcal/mol) showed higher binding scores than the reference drug pentamidine ( − 8.10 kcal/mol), while N3 and N7 – N12 showed higher predicted pIC 50 than the templates. Also, the pharmacokinetics properties prediction revealed that the newly designed compounds, obeyed the Lipinski’s rule for oral bio-availability, showed high human intestinal absorption (HIA), low synthetic accessibility score, CNS and BBB permeability and were pharmacologically active. Conclusions: The activities of the various maleimides were predicted excellently by the built QSAR model. Based on the pharmacokinetics and molecular docking studies therefore, the newly designed compounds are suggested for further practical evaluation and/or validation as potential drug candidates for the treatment of leishmaniasis.


Background
Leishmaniasis is among those diseases which are classified as neglected tropical diseases. It was reported as the second most prevalent parasite infection after malaria (Tonelli et al. 2018). Leishmaniasis is known to affect a population of about 350 million people majorly living in the tropical impoverished regions of the world, with approximately 12 million people been infected annually (Baquedano et al. 2016). It is caused by a group of protozoan parasites of the genus leishmania, which are responsible for the three common clinical types: cutaneous, mucocutaneous and visceral leishmaniasis (VL) among which VL is the most fatal if left untreated (Keurulainen et al. 2018). VL is majorly caused by Leishmania Donovani (L. donovani) and L. infantum, inherited through the bite of the female sand fly (Abdelhameed 2018). The challenge remains that only few drugs are available for the treatment of leishmaniasis including Pentostam (contains heavy metal antimony), pentamidine, amphotericin B, paromomycin and miltefosine with varying shortcomings such as high cost, toxicity and less efficacy among others (Fan et al. 2018). For example, fever, chill, nephrotoxicity, hypokalemia and myocarditis are the major adverse effects associated with amphotericin B (Ghorbani and Farhoudi 2018). Miltefosine has been linked to contraindication in pregnancy and mandatory contraception for women in child bearing age (Ghorbani and Farhoudi 2018). The major safety concern associated with pentamidine is induction of insulin-dependent diabetes mellitus, while ototoxicity and nephrotoxicity are associated with paromomycin (Seifert 2011). Also, the parasites which are increasingly becoming drug-resistant pose a great threat to the treatment of leishmanial infections (Fan et al. 2018), and more so that it has received only a little global attention compared to other infections such as malaria, tuberculosis and cancer. Therefore, the development of new medicines with clinical attributes that overcome all of these drawbacks is highly necessary.
In this study, a computer-aided drug design approach was utilized to propose some set of new compounds for further evaluation as leishmania inhibitors. In the process of new drug development, computational tools play a critical role because it saves time and cost and has proven to be more effective than the crude traditional methods (Lawal et al. 2021). Notable computational tools include quantitative structure activity relationship (QSAR), molecular docking, molecular dynamics, pharmacokinetics and homology modeling among others (Adeniji et al. 2019;Ibrahim et al. 2020Ibrahim et al. , 2021. The knowledge of QSAR helps in establishing a relationship between various molecular structures of molecules and their experimental activities (Adeniji et al. 2019). Molecular docking simulation is a computer-aided virtual screening method which probes the binding of ligands in the active sites of the protein target using a valid docking tool . Pharmacokinetics analysis is important in the preclinical study of new drug compounds in order to ascertain how such drug compounds affect the living organism when administered. Some of the most important pharmacokinetics properties to be determined during preclinical testing include absorption, distribution, metabolism, excretion and toxicity (ADMET) (Lawal et al. 2021;Ibrahim et al. 2021). Physicochemical properties such as molecular weight, topological polar surface area (TPSA), lipophilicity, water solubility, hydrogen bond donors, and hydrogen bond acceptors are necessary to predict drug's likelihood of being orally bioavailable (Lipinski et al. 2001). The choices of molecules for oral bioavailability have been guided by several rules such as the Lipinski's 'rule of 5' (RO5), Veber rule, Ghose rule, Egan and Muegge (Sun et al. 2020).
Some therapeutic protein targets of Leishmania spp available in the protein data bank for docking with drug ligands include Ribokinase from L. donovani (5ZWY) (. Gatreddi et al. 2019), L. major mitochondrial ribosome (7ANE) (Soufari et al. 2020), pyridoxal kinase from L. donovani in complex with ADP and pyridoxine (6K91) ) and L. donovani UMP synthase (3QW4) (French et al. 2011). Pyridoxal kinase (PdxK) from L. donovani is an interesting drug target; it was reported to catalyze the phosphorylation of the 5′ hydroxyl group of pyridoxal to form pyridoxal-5′-phosphate, the biologically active form of vitamin B6 . PdxK was in a recent study shown to be significant for parasite growth and key to host infection (Kumar et al. 2018). Known anti-malaria medicines like chloroquine and primaquine had earlier shown inhibition against PdxK (Kimura et al. 2014). Therefore, PdxK represents a promising target for designing novel anti-leishmanial agents. Structure-based drug design (SBDD) otherwise referred to as the direct drug design is a critical part of industrial drug discovery processes, and a useful tool in academic researches. This is because it deals with the knowledge of the 3D structures of the target enzymes and the related diseases at the molecular level (Batool et al. 2019). The method comprises of steps such as obtaining useful information on the 3D structure of the target protein determined via X-ray crystallography; NMR spectroscopy or homology modeling; retrieval and preparation of such enzyme target; and the actual design of novel compounds (Abdullahi et al. 2022).
Maleimides have been reported as antimicrobial compounds (Chen et al. 2015;Li et al. 2012;Shen et al. 2013) and also have shown excellent enzyme inhibitory activities (Silvia and Maria 2005;Slavica et al. 2007). Additionally, maleimides had earlier been described as anti-inflammatory (Nara et al. 2010), anti-cancer (Khan et al. 2004) and anti-anxiety (Jerzy 2003). In order to exploit the anti-leishmanial effect of maleimides therefore, Fan et al. (2018) synthesized a series of maleimides and reported their biological activities against L. donovani. Consequently, this study was undertaken to develop a QSAR model for predicting with precision the activities of some maleimides as potent anti-leishmanial agents, conduct a virtual docking screening with different protein receptors, carry out a structure-based design of new maleimides and subjecting same to pharmacokinetics and molecular docking studies in order to evaluate their drug-likeness properties and binding interaction pattern, respectively.

Data collection
A series of 28 maleimides with reported anti-leishmanial activities (IC 50 in µg/mL) were obtained from the literature (Fan et al. 2018). The bioactivities (IC 50 ) were first converted from g/L to Molar (M) using Eq. (1) and thereafter to logarithmic scale (pIC 50 ) using Eq. (2) (Wang et al. 2020). The molecular structures of the various maleimides with their observed IC 50 and pIC 50 values are shown in Table 1.

Molecular geometry optimization
The two-dimensional (2-D) structures of the various compounds were drawn using ChemDraw Ultra (version 12.0.2) and saved as MDL mol file format. The resulting files were then fed separately onto the Spartan'14 (version 1.1.4) graphical user interface in threedimensional (3-D) structural form. The 3-D structures were first subjected to energy minimization and molecular mechanics force field (MMFF) in order to minimize their chemical structures and to remove tension energy of the molecules' conformation. The main optimization process was then conducted using density functional theory (DFT) with Becke's three-parameter read-Yang-Parr hybrid (B3LYP) option and utilizing the 6-31G basis set. The thoroughly optimized structures were then saved in SD file and PDB formats for use in descriptor calculation and molecular docking, respectively (Wang et al. 2020;Li et al. 2004).

Generation of molecular descriptors
The resulting data in SD file format obtained earlier from the optimization process were imported into the Pharmaceutical Data Exploration Laboratory (PaDEL)-Descriptor software (version 2.20) in order to calculate the molecular descriptors for all twenty-eight (28) compounds (Lawal et al. 2021).

Dataset pretreatment and division into training and test sets
The Drug Theoretical and Chem-informatics Laboratory (DTC Lab)-based software GUI 1.2 was used to remove non-informative descriptors from the generated descriptor pool (Adeniji et al. 2020). The pretreated data were then divided into the modeling train set data and external evaluation test set data in the ratio of 70:30, respectively, with the help of DTC Lab-derived software which utilizes the Kennard Stone method for dataset division (Kennard and Stone 1969). The splitting of dataset into training and test sets was based on the closeness of the representative points of the test set to the representative points of the training set in the multidimensional descriptor space (Ugbe et al. 2021).

Model building
The genetic function approximation (GFA) as a statistical technique in the Material Studio software (version 8.0) was used to generate the models based on multi-linear regression (MLR) approach. GFA was used to obtain the optimum descriptor combinations constituting the QSAR models, while MLR helps to establish the relationship between the biological activities, pIC 50 (dependent variable) and the molecular descriptors (independent variables) . The multi-linear regression equation assumes the following form (Eq. 3) (Adawara et al. 2020): where Y represents the dependent variable; 'k's and 'x's represent, respectively, regression coefficients and independent variables; 'C' equals intercept or regression constant.

Assessment of model quality
Internal validation assessment of the built model was carried out on Material studio by GFA approach using the Friedman formula, correlation coefficient (R 2 ) and crossvalidation coefficient (Q 2 cv), and also by the Y-randomization test in order to show how well the model predict the activity values of the model building compounds (Adeniji et al. 2018;Adawara et al. 2020). The model's predictive power was also assessed externally to show if the model could predict the activity values of the test  (Isyaku et al. 2020). Furthermore, the data were subjected to the Golbraikh and Tropsha acceptable model criteria using the MLR-plusValidation tool (version 1.3) (Roy et al. 2013;Edache et al. 2020). The various equations and parameters used for the model validation are presented in Table 2.

Evaluation of descriptors' mean effects
The mean effect (ME) value shows the relative contribution of each descriptor in a model, defined as (Eq. 10): (10) ME = B j n i D j m j B j n i D j   where β j is the coefficient of the descriptor j in the model, D j is the value of each descriptor in the data matrix for each molecule in the training set, m is the number of the descriptor that appears in the model, n is the number of molecules in the training set (Abdullahi et al. 2019).

Variance inflation factor
The degree of multi-co-linearity or correspondence between the descriptors is measured by the variance inflation factor (VIF), usually defined as (Eq. 11): where R 2 is the correlation coefficient of the multiple regressions between the variables within the model. VIF values of 1 indicate no inter-correlation exists for each variable; for VIF in the range of 1-5, the related model is acceptable; and if VIF is greater than 10, the related model is unstable and unacceptable (Abdullahi et al. 2019).

Evaluation of the model's applicability domain
Evaluating the applicability domain (AD) of a QSAR model is important to ascertain the reliability and robustness of the built QSAR model. AD provides one the chance to estimate the uncertainty in the prediction of compounds based on their similarity with the training set compounds, used in the model building (Tropsha et al. 2003). The leverage approach was used to describe the AD of the developed model. The leverage (h) of a particular chemical compound is defined thus (Eq. 12): where X = m × k descriptor matrix of the training set compound, X T = transpose matrix of X.
The warning leverage (h*) which is the range of values used to check for influential molecule or outlier is defined below (Eq. 13): in the training set, Y training = mean experimental activity of the training set, Y exp experimental activity in the training set, Y pred predicted activity in the training set, n number of compounds in the training set, cR 2 p Y-randomization coefficient, R correlation coefficient for Y-Randomization, R r average 'R' of random models, Ypred test predicted activity of test set, Yexp test experimental activity of test set, r 2 square correlation coefficients of the plot of experimental activity versus predictedactivity values, r o 2 square correlation coefficients of the plot of experimental activity versus predicted activity values at zero intercept, rʹ o 2 square correlation coefficients of the plot of predicted activity versus experimental activity at zero intercept, kʹ slope of the plot of predicted activity against experimental activity at zero intercept

Parameter Equation Eq Significance Threshold value
Internal validation Allows for the best fitness score to be obtained Measures the degree of fitness of the regression equation Ensures model's stability and reliability ≥ 0.5 Indicates a high internal predictive power ≥ 0.5 The coefficient of determination ( cR 2 p ) of Y-Randomization This is for a confirmation that the QSAR model built is strong and not created by chance Measures the ability of the model to predict activity values of external set of compounds ≥ 0.6 Golbraikh and Tropsha acceptable model criteria Assess the robustness and stability of the model where m number of training set compounds, j number of descriptors in the model. A plot of the standardized residuals against leverages otherwise called the William's plot was used to evaluate the significant area in the model's chemical space. As a rule, compounds which fall within this area on the plot are the approved predicted compounds (Adeniji et al. 2020;Veerasamy et al. 2011).

Virtual docking screening
Molecular docking investigation was conducted separately between five (5) different protein receptors of L. donovani and all 28 maleimides using the Auto Dock Vina of PyRx v software tool. To revalidate the docking results, the reference drug (pentamidine) was equally docked onto the same binding pockets of the 5 receptors. This screening was carried out in order to identify the best therapeutic protein target for the maleimide molecules. The various receptors in 3-D form were retrieved from the protein data bank in the PDB file format and fully prepared on the Biovia Discovery Studio Visualizer by excluding water molecules and co-crystallized ligand enclosed within the protein structure. The 3D structures of all the ligands (maleimides) were saved in PDB file format after the optimization process using Spartan 14 (Adeniji et al. 2019;Adawara et al. 2020). The various receptors used in the virtual docking screening are described in Table 3.

Structured-based drug design
In this study, twelve (12) compounds (N1-N12) were designed via SBDD using the template compounds (14, 21 and 24) basically by addition of varying substituent groups at the meta, para and ortho positions of the phenyl ring system in the various template structures. N1-N3 were designed out of compound 14, N4-N6 out of 21 and N7-N12 via compound 24. The molecular structures of the newly designed compounds were drawn using the ChemDraw Ultra and subjected to molecular geometry optimization according to the procedures earlier reported herein for molecular geometry optimization. The optimized structures were then saved in SD and PDB files format for further analyses involving molecular descriptor calculations and molecular docking study (Abdullahi et al. 2022). Standard docking was performed between the newly designed maleimides and the preferred target receptor (PdxK), using the Auto Dock Vina of PyRx v software tool, while the Biovia Discovery Studio tool was used to visualize the resulting pharmacological interactions as earlier reported herein for virtual docking screening.

Prediction of pharmacokinetic properties
Pharmacokinetics properties prediction constitute an absolutely necessary stage in drug discovery's early phase because only molecules with good drug-likeness properties and excellent ADMET profiles advance into the preclinical research phase (Lawal et al. 2021). Therefore, all the newly designed compounds were investigated for their drug-likeness and ADMET properties using the online web servers; http:// www. swiss adme. ch/ index. php and http:// biosig. unime lb. edu. au/ pkcsm, respectively. The Lipinski's RO5 is a widely used criterion for oral bioavailability. Hence, the tested compounds would be assessed for oral bioavailability using the RO5 criteria (Lipinski et al. 2001).

QSAR study
The results of theoretical study (QSAR) conducted on the twenty-eight (28) maleimides are presented in Figs. 1, 2, 3 and Tables 4, 5, 6, 7, 8. The built QSAR model (Eq. 14) is composed of four (4) descriptors described clearly in Table 4. In order to ascertain the stability, robustness, reliability and predictive power of the built QSAR model, it was subjected to both internal and external validation tests, and the results are presented in Table 5. The calculated descriptors, experimental activities (pIC 50 ) of the various compounds, together with those predicted by the model as well as their residuals, are presented in Table 6. Also, Fig. 1 shows the plot of predicted activity values against those of experimental activity for both training and test sets. A further plot of standardized residuals against experimental activities was obtained and is presented in Fig. 2.     Table 7. Another significant validation test is the Y-randomization test, which was also performed in order to confirm that the QSAR model is strong and not created by chance. The resulting parameters are presented in Table 8. Also, the scatter plot of the standardized residuals versus the leverages (William's Plot) obtained to  ascertain the model's applicability domain is shown in Fig. 3.

Virtual docking screening
The results (binding affinities) of the docking simulation conducted between five (5) protein receptors of L. donovani and the studied maleimides, as well as the reference drug (pentamidine) are reported in Table 9. Also, the prepared 3-D structure of PdxK is shown in Fig. 4, while the pharmacological interactions between the receptor (PdxK)'s amino acid residues and the selected templates (14, 21 and 24), as well as the reference drug (pentamidine) as compiled through Discovery Studio Visualizer, are presented in Table 10. Additionally, the 2D and 3D views of the binding interactions of compounds 14, 21, 24, and pentamidine with the binding site of PdxK receptor as adapted from the Discovery Studio Visualizer, are shown in Figs. 5, 6, 7 and 8, respectively.

Structured-based drug design
The 2D molecular structures, predicted pIC 50 and binding energies of the newly designed compounds as well as the reference drug (pentamidine) are reported in Table 11. Also, the binding interaction profiles of the target receptor (PdxK) with the newly designed compounds are presented in Table 12, while the 2D and 3D views of the binding interactions of N7 and N12 with the receptor's amino acid residues are shown in Figs. 9 and 10, respectively.

Pharmacokinetics properties prediction
Results of the pharmacokinetics investigation conducted on the newly designed compounds are presented in Tables 13, 14, while Fig. 11 shows their Boiled Egg's representation.

QSAR study
A combined GFA and MLR approach led to the selection of four (4) descriptors in the generation of the QSAR model. The built model (Eq. 14) was found to excellently satisfy the requirement for a reliable QSAR model. The low residual values between the experimental and predicted activities as shown in Table 6 indicate a high predictive strength for the QSAR model. The R 2 values of 0.801 and 0.892 for training set and test set, respectively, as obtained from the plot of exp. pIC 50 against pred. pIC 50 in Fig. 1 compare perfectly well with those obtained from GFA (0.8013 and 0.892) and MLRplusValidation analysis (0.8013 and 0.892) as reported in Table 5.  The close grouping of points along the line of best fit in Fig. 1 shows a very good correlation between the experimental and predicted activity values, indicating that the built model is reliable and robust. The random spread of standardized residuals on both sides of the zero mark in Fig. 2 is an indication that the built model is free of any systematic error.
Furthermore, the low correlation coefficients (less than 0.50) which exist between each pair of the descriptor in the built model (Table 7) indicate no inter-correlation between each descriptor. Similar result was also obtained elsewhere by Adeniji et al. (2018) and Abdullahi et al. (2022). The VIF values ranging from 1 to 5 for all 4 descriptors as reported in Table 7 showed that  the descriptors were statistical orthogonal and the built model was statistically substantial, an indication of the stability and acceptability of the built QSAR model. Similar observation was reported by Adeniji et al. (2019). The values of the absolute t-statistics greater than 2 for each descriptor show that the selected descriptors were good (Adeniji et al. 2018). Also, the evaluated p-values for the various descriptors in the model at 95% confidence level were less than 0.05 as shown in Table 7. Therefore, the alternative hypothesis which asserts that a relationship exists between the descriptors used in generating the model and the compounds' inhibitory activities at p ˂ 0.05 holds. Additionally, the values of the mean effect (ME) reported in Table 7 provide vital information on the effect and degree of each descriptor's contributions in the built model. The magnitudes and signs of ME values signify their individual strength and direction on the molecules' inhibitory activities. All the descriptors except ATSC6v have positive ME, meaning that increase or decrease in their values will lead to an increase or decrease in the anti-proliferative activities, respectively. Increasing the values of ATSC6v will lead to a decrease in the inhibitory activities because of its negative ME value. GATS5c with the highest ME value has the greatest influence on the molecules' inhibitory activities. GATS5c is Geary autocorrelation of lag 5 weighted by gasteiger charge, which has a positive ME suggested to contribute positively to anti-leishmanial activity. The gasteiger charge is a physicochemical property calculated for every atom in the molecule and is the Geary coefficient (Mahmud et al. 2020).
The low values of R 2 and Q 2 obtained from the random reshuffling (Table 8) inferred that the built model is stable, robust and reliable. The value of coefficient for Y-randomization, cR 2 p (0.664372) greater than 0.50, supports the claim that the built model is powerful and not inferred by chance. The William's plot (Fig. 3) clearly shows that all the compounds fall within the square area ± 2.5 of standardized cross-validated residual. It can therefore be inferred that no outlier is present in the dataset. However, compound 1 was found with leverage value greater than the calculated warning leverage (h* = 0.75) and was said to be an influential molecule.

Virtual docking screening
Binding energies of the protein-ligand (drug) interactions are important to describe how well the drug binds to the target macromolecule. The negative value of the binding energy change shows the spontaneity of the binding process and how well ligands can fit into the target protein pocket to form the most energetically stable drug receptor (Ugbe et al. 2021). Among the studied receptors, pyridoxal kinase (PdxK) receptor (pdb id: 6K91) had relatively shown the strongest interaction with the various compounds as shown by the higher binding energy values associated with this receptor (Table 9). Consequently, PdxK was selected as the target receptor of interest in this study and subsequent discussions pertaining protein-ligand interactions shall be based on it. Also, among the 28 maleimides studied, compounds 14, 21 and 24 bind more strongly to PdxK with the highest  reported binding energies of − 7.7, − 7.7 and − 7.8 kcal/ mol, respectively. These compounds were equally well predicted by the built QSAR model with low residual values and contained within the model's applicability domain (Table 6 and Fig. 3). Therefore, 14, 21 and 24 were selected as template molecules for designing new compounds with improved binding scores and pharmacological properties. As seen from Table 10 and Figs. 5, 6, 7, 8, compound 14 was observed to have interacted well with the binding site of the PdxK receptor through three (3) conventional hydrogen bonds, one (1) carbon-hydrogen bond, one (1) π-anion, one (1) π-π stacked and one (1) pi-alkyl interactions. One of the carbonyl oxygen atoms on the pyrrole ring system formed 3 conventional hydrogen bonds with 2 THR-229 at distances of 2.41 Å and 2.48 Å, and GLY-228 at a distance of 2.37 Å. It also formed a carbonhydrogen bond with ASP-231 at a distance of 3.44 Å. Others include π-anion between the pyrrole ring system and ASP-124, pi-pi stacking between the phenyl  ring and aromatic ring of TYR-85, and π-alky interaction between the phenyl ring and VAL-19. Similarly, Compound 21 binds well into the binding pockets of the PdxK receptor via five (5) conventional hydrogen bonds, a carbon-hydrogen bond, one π-anion, two π-alkyl, and one alky interactions. One of the carbonyl groups oxygen ASN asparagine, ASP aspartic acid, GLN glutamine, GLU glutamic acid, GLY glycine, HIS histidine, ILE isoleucine, LEU leucine, LYS lysine, MET methionine, PHE phenylalanine, SER serine, THR threonine, TYR tyrosine, VAL valine The ligand also formed a carbon-hydrogen bond with THR-227 at a distance of 2.67 Å. Others include alkyl interaction with VAL-121, π-anion with ASP-124, and π-alkyl with LYS-187 and LEU-257. Distinguishably, the binding interactions of Compound 24 with the receptor were characterized only by electrostatic and hydrophobic interactions, and without hydrogen bond interactions. This may be attributed to the steric hindrances posed by the side chain methyl groups on the phenyl ring and the chloro groups on the pyrrole ring which shields the carbonyl groups from interacting with amino acid residues to form hydrogen bonds. The observed interactions were dominated by the hydrophobic interactions (alkyl type) with the exception of ILE-261 which binds by electrostatic interaction to the π electron system of the benzene ring of compound 24 via its alkyl group. Lastly, the reference drug (pentamidine) was equally docked onto the  binding pocket of the pyridoxal kinase receptor in order to provide insight into their binding interaction mode and for validation purpose. Pentamidine interacted with PdxK via two (2) conventional hydrogen bonds with ASN-151 and SER-12 at a distance of 1.95 Å and 2.28 Å, respectively. It equally formed a carbon-hydrogen bond with ASP-231 at a distance of 3.40 Å. Others are electrostatic and hydrophobic interactions including π-alkyl interactions with VAL-19 and LYS-187, π-cation with TYR-85 and π-π stacked with TYR-85.

Structured-based drug design
As reported in Table 11, the predicted pIC 50 values of most of the designed compounds (N3, N7-N12) were greater than those of their corresponding template molecules, an evidence of being more biologically active molecules than their templates. Also, all the newly designed compounds showed higher docking scores compared to the template molecules, while only N7 and N12 with binding energy values of − 8.9 and − 8.5 kcal/mol showed better docking scores than that of the reference drug (pentamidine). This is a clear indication of how well these new maleimides would interact with the target enzyme. Because N7 and N12 bind more strongly to PdxK than pentamidine, their interactions with the binding pocket of the protein target shall be discussed further. N7 was observed from Fig. 9 to bind excellently with the active site of PdxK via six (6) conventional hydrogen bonds, one (1) π-alkyl hydrophobic interaction and one (1) π-anion electrostatic interaction. The carbonyl groups oxygen atoms play a vital role in hydrogen bond formation just as observed with the template molecules. One of the carbonyl groups formed hydrogen bond with ASN-151 and LYS-187 via its oxygen at distance of 2.32 Å and 2.11 Å, respectively. The other carbonyl group formed a hydrogen bond with GLY-228 at a distance of 2.40 Å and 2 hydrogen bonds with THR-229 at a distance of 2.38 Å and 2.40 Å. It also formed an additional hydrogen bond with GLN-258 at a distance of 1.98 Å via a hydroxyl group at the para-position of the outer phenyl ring system. Others include π-alkyl hydrophobic interaction with LEU-198 via the compound's π electrons of the outer benzene ring, and π-anion electrostatic interaction with ASP-124 via the compound's π electrons of the pyrrole ring system. The interaction of N12 with the protein target as observed from Fig. 10 was via three (3) conventional hydrogen bonds, three (3) carbon-hydrogen bonds and one (1) π-donor hydrogen bond. The observed hydrophobic interactions include one each of alkyl and π-alkyl interactions, while the electrostatic interactions include one each of π-anion and π-cation interactions. Also, observed was a halogen interaction with the receptor. The observed conventional hydrogen bonds were formed by the interaction of its carbonyl groups oxygen with LYS-187 and THR-229 at distance of 2.06 Å and 2.46 Å, respectively, and between its nitro group oxygen and SER-188 at a distance of 2.64 Å. The carbon-hydrogen bonds were formed by its trifluoromethyl group's interaction with GLY-230, LYS-187 and THR-229 at interaction distance of 2.71 Å, 3.04 Å and 2.38 Å, respectively. Also observed was a π-donor hydrogen bond between the π-electron of the pyrrole ring system and ASN-151 at a distance of 3.19 Å. The observed hydrophobic interactions include alkyl interaction with LYS-187 and VAL-121, and π-alkyl interaction with LYS-187. Furthermore, the observed electrostatic interactions comprises of π-cation between LYS-187 and its benzene ring π-system, and π-anion between ASP-124 and its pyrrole ring π-system. Therefore, these compounds have demonstrated the potentials to arrest the target receptor (PdxK) an absolutely necessary factor whose inactiveness is dangerous to the viability of the parasite (L. donovani).

Pharmacokinetics properties prediction
According to the Lipinski's rule for oral bioavailability, a drug molecule is more likely to have poor absorption or permeation when it has hydrogen bond donors (HBD) of greater than 5, hydrogen bond acceptors (HBA) > 10, molecular weight (MW) > 500 and lipophilicity (MLOGP > 4.15 or WLOGP > 5) (Lipinski et al. 2001). Molecules that satisfy at least three out of the four requirements are said to obey the Lipinski's rule for oralbioavailability (Lawal et al. 2021).
As seen from Table 13, all the designed molecules perfectly obeyed the Lipinski's rule by showing no violation. Also, the reported values of topological polar surface area (TPSA) for all molecules were less than 140 Å 2 , beyond which a molecule may exhibit poor gastrointestinal absorption. Additionally, the synthetic accessibility scores of all tested molecules are in the easy portion (˂ 5.00), indicating their easy laboratory synthesis. The estimated water solubility (Log S) ranges from moderately soluble (N5, N8 and N12) to soluble (the remaining 9 compounds). The Boiled Egg representation in Fig. 11 provides for an intuitive evaluation of passive gastrointestinal absorption and brain barrier permeability as a function of the position of the molecules in the WLOGP against TPSA plot (Daina et al. 2017). As seen from Fig. 11, all the designed molecules were represented in red dots, indicating that they were predicted not to be effluated from the central nervous system by P-glycoprotein. P-glycoprotein is an enzyme which acts as a biological barrier by extruding toxins and xenobiotics, including drugs out of cells. Also, six (6) molecules (N1, N6, N7, N9, N10 and N12) were located in the Boiled Egg's yolk, meaning that these molecules were predicted to passively permeate through the blood-brain barrier (BBB), while the remaining six (6) molecules were located in the Boiled Egg's white which is an indication that they were predicted to be passively absorbed by the gastrointestinal tract.
The predicted ADMET properties in Table 14 showed that the human intestinal absorption (HIA) was high (> 70%) for all newly designed compounds. Drug molecules are said to be poorly distributed to the brain through the blood-brain barrier (BBB) and considered as unable to penetrate the central nervous system (CNS), when the values of the logarithmic ratio of brain to plasma drug concentration (logBB) are less than − 1, and the blood-brain permeability-surface area product (logPS) are less than − 3, respectively. Consequently, all the newly designed compounds were predicted to cross the BBB except N7, N10 and N12 with predicted logBB values of less than − 1 (Table 14). Also as predicted in Table 14, all the molecules showed CNS permeability, i.e., logPS > − 3. Furthermore, all the studied molecules except N7 and N11 are substrates of the Cytochrome P450 enzyme (CYP-3A4), an important enzyme for drug metabolism in the body, with none of the compounds inhibiting the enzyme. The degree of drug elimination from the body is measured by the drug's total clearance, which is within the accepted range for these newly designed compounds. Additionally, the predicted values of maximum recommended tolerated dose (MRTD) for all molecules are included in Table 14. MRTD value of less than or equal to 0.477 log (mg/kg/day) is considered low, and high if greater than 0.477 log (mg/kg/day). The overall predicted drug-likeness and ADMET properties put these molecules on an excellent pharmacokinetics profile, and more so that they are orally bio-available. Therefore, the newly designed molecules are suggested for practical evaluation and/or validation in the laboratory as potential drug candidates for the treatment of leishmaniasis.

Conclusions
In this study, a four-descriptor QSAR model was developed with a series of twenty-eight (28) maleimides, which was used to excellently predict their antileishmanial activities, and those of the newly designed compounds. The virtual docking screening conducted between the 28 maleimides and five (5) target enzymes, identified pyridoxal kinase (PdxK) receptor as a biological target of interest. Compounds 14, 21 and 24 were used as templates for designing twelve (12) new maleimides, because of their relatively higher binding energies of interaction with PdxK. The newly designed compounds interacted more strongly with the target site than the chosen template molecules, but only N7 and N12 were better than the reference drug (pentamidine) in this regard. The predicted pharmacological interaction profiles of the designed compounds generally showed a good fitting into the target site cavity. Also, all the newly designed compounds were said to be orally bio-available, and showed good drug-likeness and ADMET properties. Hence, these new molecules have demonstrated the potential to arrest PdxK enzyme, which could inhibit the growth and viability of the parasite (L. donovani). It is therefore recommended that, these new compounds be developed and further evaluated as potential drug candidates for the treatment of leishmaniasis.