A combined 2-D and 3-D QSAR modeling, molecular docking study, design, and pharmacokinetic profiling of some arylimidamide-azole hybrids as superior L. donovani inhibitors

Leishmaniasis is one of the neglected tropical diseases which is prevalent in the tropical regions of the world most especially in Africa. It is caused by the Leishmania species and transmitted to humans majorly through the bite of the female sandfly. This study was carried out in support of the continuous search for new drug molecules effective enough for the treatment of leishmaniasis, and which have very limited side effects. This study was focused on a combined 2-D and 3-D QSAR modeling of thirty-six arylimidamide-azole hybrids, molecular docking study, design, and pharmacokinetic analysis of some selected and newly designed arylimidamide-azole analogs. The density functional theory (DFT) with B3LYP and 6-31G** basis set was employed for the geometry optimization of the various compounds. Genetic function approximation (GFA) and multi-linear regression (MLR) approaches were used for the 2-D QSAR model building, while the fractional factorial design (FFD) and uninformative variable elimination-partial least square (UVEPLS) were employed for building the 3-D QSAR model. Pyridoxal kinase (PdxK) receptor (PDB: 6K91) was the target protein of interest in this study. The built 2-D and 3-D QSAR models were found to satisfy the requirement of both internal and external validation tests as follows: 2-D QSAR; R2 = 0.9614, R2adj = 0.9513, Q2cv = 0.9350, R2test = 0.6766 and cRp2 = 0.8779, and for 3-D QSAR (UVEPLS at PC = 5); R2 = 0.9839, Q2LOO = 0.7539 and Q2LTO = 0.7367. The CoMFA steric and electrostatic field contributions were 68.96% and 31.04%, respectively. All the designed analogs showed higher predicted activities than that of the template (36). Also, the new compounds showed higher binding affinities (MolDock scores) than that of the reference drug pentamidine (− 141.793 kcal/mol), with 36e showing the highest negative MolDock score of − 208.595 kcal/mol. Additionally, these newly designed compounds were said to be orally bioavailable (excluding 36f and 36g that violated 2 of the Lipinski’s provisions), with perfect intestinal absorption, less difficult to synthesize, AMES toxicity free, and showed strong interactions with the target. The newly designed compounds especially 36e have shown a marked pharmacological improvement over the template molecule and are therefore recommended for further practical evaluation as superior pyridoxal kinase inhibitors.

Page 2 of 24 Ugbe et al. Bulletin of the National Research Centre (2022) 46:189 Background Leishmania species is a group of protozoan parasites responsible for leishmaniasis, a neglected tropical disease known to affect a global population of up to 350 million people majorly in the tropical regions (Al-Tamimia et al. 2019;Baquedano et al. 2016). Three common clinical manifestations were earlier reported as cutaneous, mucocutaneous, and visceral leishmaniasis (VL) (Díaz et al. 2019). Among these forms, cutaneous is the most common, while VL is the most fatal (Keurulainen et al. 2018). Specific species causing VL include Leishmania Donovani (L. donovani) and L. infantum, transmitted by the insect vector (female sandfly) (Ugbe et al. 2022a, b). Only a limited number of drugs are readily available for the treatment of leishmaniasis such as pentamidine, amphotericin B, pentostam, paromomycin, and miltefosine, none of which is reported free from associated drug adverse effects, low efficacies, high cost, etc. (Fan et al. 2018). For example, amphotericin B has been linked to fever, nephrotoxicity, myocarditis, chill, and hypokalemia (Ghorbani and Farhoudi 2018). Induction of insulin-dependent diabetes mellitus has been the major drawback of pentamidine, while paromomycin is not free from nephrotoxicity and ototoxicity (Seifert 2011). Also, most organisms had in recent times shown constant resistance to some of these medicines, a case which has greatly threatened the effective treatment of leishmaniasis (Fan et al. 2018). Additionally, the disease has received less attention compared to other infections such as hepatitis, malaria, AIDS, diabetes, tuberculosis, and cancer. As such, it has become imperative to develop new medicines with attributes that overcome all of these shortcomings.
In the continuous search for new effective drug compounds, in silico approaches like the molecular docking studies, quantitative structure-activity relationship (QSAR), homology modeling, molecular dynamics, and pharmacokinetics analysis among others play a very crucial role because it is cost-effective, save time and proves effective than the crude methods (Adeniji et al. 2019;Ugbe et al. 2021). The connection between a compound's structural properties (independent variables) and its biological activities (dependent variable) is usually established by QSAR (Adeniji et al. 2019). Molecular docking on the other hand helps to investigate the binding of ligands with the active sites of the target receptor . The pharmacokinetic study is a crucial step in the preclinical phase of new drug molecules, which helps to ascertain how such drug compounds affect the living organism when administered (Lawal et al. 2021;Ibrahim et al. 2021). The choice of drug molecules for oral bioavailability has been widely guided by the Lipinski rule of five (RO5) or the Pfizer rule which considers physicochemical properties such as molecular weight (MW), lipophilicity, hydrogen bond donors (HBDs), and hydrogen bond acceptors (HBAs) as necessary factors in predicting drug's likelihood of being orally bioavailable (Lipinski et al. 2001). Also important for evaluation as part of a pharmacokinetics study are absorption, distribution, metabolism, excretion, and toxicity (ADMET) (Ugbe et al. 2022a, b).
Pyridoxal kinase (PdxK) from L. donovani (PDB: 6K91) is an interesting enzyme previously reported to catalyze the phosphorylation of the 5′ hydroxyl group of pyridoxal to form pyridoxal-5′-phosphate, an active form of vitamin B6 . PdxK is an important enzyme for parasite growth and key to host infection (Kumar et al. 2018). Chloroquine and primaquine which are well-known anti-malaria drugs were earlier reported to inhibit PdxK (Kimura et al. 2018). Therefore, PdxK serves as a promising protein target for designing new leishmanial inhibitors. Arylimidamideazole compounds incorporate both arylimidamide and azole moiety in their structures, which were reported to possess the advantage of interacting with the targets of both classes of compounds and which also show improvement in the in vivo pharmacokinetics and/ or pharmacodynamics of these classes of molecules (Abdelhameed et al. 2020).
Therefore, this study was focused on a combined 2-D and 3-D QSAR modeling for predicting the activities of some arylimidamide-azole hybrids, performing molecular docking studies, and design of new analogs, while subjecting same to pharmacokinetics analysis to evaluate their drug-likeness and ADMET properties.

Data acquisition
A series of thirty-six arylimidamide-azole hybrids with reported biological activities (IC 50 in µM) against L. donovani intracellular amastigotes were sourced from the literature (Abdelhameed 2018). The various bioactivity (IC 50 ) values were separately converted to pIC 50 using Eq. (1) (Ugbe et al. 2022a, b). The molecular structures of the various arylimidamide-azole hybrids Keywords: Arylimidamide-azole, DFT, Leishmaniasis, Molecular docking, Pharmacokinetics, 2-D QSAR, 3-D QSAR    with their observed anti-leishmanial activities, including the molecular structure of the reference drug (pentamidine), are shown in Table 1.

Molecular geometry minimization
The molecular structures of all the compounds were drawn using the ChemDraw Ultra, saved as MDL molfile format, and thereafter imported separately onto the Spartan'14 Graphical User Interface while enabling the auto conversion of 2-D models to 3-D. The imported molecules were initially subjected to energy minimization and then saved in Spartan file format. The resulting structures were then fully optimized first by using molecular mechanics force field (MMFF) and thereafter density functional theory (DFT) with Becke's three-parameter read-Yang-Parr hybrid (B3LYP) option and utilizing the 6-31G** basis set. The optimized structures were then saved as SD files and PDB formats for subsequent use in descriptor generation and docking, respectively Ugbe et al. 2021).

2-D QSAR model building
The 2-D QSAR model building began with the generation of molecular descriptors for all thirty-six compounds by feeding the various molecules in SD file format into the PaDEL-Descriptor software (Lawal et al. 2021). Data pretreatment using the DTC-Lab pretreatment tool (GUI 1.2) was carried out to remove uninformative descriptors from the generated descriptors pool (Adeniji et al. 2020). With the aid of the DTC-Lab data division tool which uses the Kennard Stone method, the resulting data were then partitioned into a training set and test set data in the 70:30 ratio, respectively (Kennard and Stone 1969). The Material studio software was used for the model building by employing the genetic function approximation (GFA) to obtain the optimum descriptor combination which is contained in the QSAR model equations built based on the multi-linear regression (MLR) approach. MLR establishes the relationship between the dependent variable (pIC 50 ) and the independent variables (molecular descriptors) . The MLR equation takes the general form (Eq. 2) (Adawara et al. 2020): where 'k's and 'x's are, respectively, regression coefficients and independent variables, Y is the dependent variable, and 'C' represents intercept or regression constant.
After the model building has been completed, the next important step was to subject the built models to an internal and external validation assessment. Some parameters computed from the Material studio are useful for the internal validation such as correlation coefficient (R 2 ), adjusted correlation coefficient (R 2 adj ), and cross-validation coefficient (Q 2 cv). The Y-randomization test was also necessary to show how well the model predicts the activities of the training set compounds (Adawara et al. 2020;Ugbe et al. 2021). To ascertain the built model's ability to predict the activities of the external test set compounds, a validation assessment was carried out externally. The predictive strength of the QSAR model depends heavily on the correlation coefficient (R 2 test) for the external test set (Isyaku et al. 2020). Furthermore, the Golbraikh and Tropsha acceptable QSAR model criteria were equally considered (Roy et al. 2013;Edache et al. 2020a, b;Ugbe et al. 2022a, b). Table 2 shows the selected equations and parameters used for the validation assessment.
Furthermore, the relative contribution of each descriptor in the built QSAR model is shown by the mean effect (ME) values, defined as (Eq. 9): where β j represents the coefficient j descriptor in the model, D j equals each descriptor's value for each molecule in the training set, m represents the number of descriptors present in the model, and n is the total number of compounds in the training set (Adeniji et al. 2019).
The variance inflation factor (VIF) defined as (Eq. 10) shows the extent of inter-correlation between the descriptors.
R 2 represents the coefficient of correlation of the multiple regressions between the variables in the model. VIF equal to 1 shows that no inter-correlation exists for each variable, and for VIF values within the range of 1-5, the related model is acceptable, and where VIF is above 10, the related model is unstable and therefore cannot be accepted (Abdullahi et al. 2019).
Also, evaluating the applicability domain (AD) of a QSAR model helps to ascertain the uncertainty in the prediction of compounds based on their similarity with the model building set. It is therefore a test of the reliability and robustness of the built QSAR model (Tropsha et al. 2003). The AD of the built model was described using the leverage approach. The leverage (h) of a given molecule is defined thus (Eq. 11): where X = m × k descriptor matrix of the training set compound, and X T = transpose matrix of X. In order to check for influential molecule or outlier, it is necessary to define the warning leverage (h*) because it provides the cutoff value (Eq. 12): where m = number of training set compounds, and j = number of descriptors in the model. The plot of standardized residuals versus leverages (William's plot) was used to evaluate the area of significance within the model's chemical space. Therefore, molecules that are found within this area on the plot are said to be the approved predicted compounds (Veerasamy et al. 2011;Adeniji et al. 2020).

Molecular alignment
The alignment of molecular structures plays a critical role in 3D-QSAR modeling (Al-Attraqchi et al. 2022) as it strongly determines the predictive accuracy and statistical quality of any given 3D-QSAR model (ElMchichi et al. 2020). Several alignment methods have been reported previously such as atom-based alignment, docking-based alignment, pharmacophore-based alignment, and co-crystallized conformer-based alignment among others (Zhang et al. 2020;Al-Attraqchi et al. 2022). In this study, the atom-based alignment was adopted using the Open3DAlign (O3A) software. The atom-based method attempts to match the atoms of the various structures to be aligned with those of the template structure, based on the atom's properties such as the partial charge. Compound 32 with the highest O3A_score was chosen as the template molecule, onto which the remaining thirty-five compounds were superimposed (aligned) and used to build the 3-D QSAR model.

Internal validation
Friedman lack of fit (LOF) (3) Allows for the best fitness score to be obtained (4) Measures the degree of fitness of the regression equation (5) Ensures the model's stability and reliability ≥ 0.5 This is for confirmation that the QSAR model built is strong and not created by chance Golbraikh and Tropsha's acceptable model criteria

External validation
Assess the robustness and stability of the model

Model development (CoMFA)
The aligned structures were used for building the 3-D QSAR model using the Open3DQSAR software (Edache et al. 2020a, b). 3-D QSAR can be comparative molecular field analysis (CoMFA) concerned with only steric and electrostatic fields' contributions, or comparative molecular similarity indices analysis (CoMSIA) which deals with the steric, electrostatic, hydrophobic, hydrogen bond donor, and hydrogen bond acceptor fields' contributions, among others (ElMchichi et al. 2020).
Herein, the CoMFA model is studied. A dataset of 36 compounds was divided into a training set and a test set of 25 and 11 molecules, respectively, i.e., percentage ratio of 70:30. The steric and electrostatic molecular interaction fields (MIFs) analysis was carried out on the aligned compounds placed within a 3-D cubic lattice of grid size 1.5 Å and a 5.0 Å out gap (Tosco and Balle 2011). Variables pretreatment was carried out as follows: energy cutoff (30.0 kJ/mol), elimination of variables having constant or near-constant values, and standard deviation cutoff (level = 2.0) (Al-Attraqchi et al. 2022). The fractional factorial design (FFDSEL) and uninformative variable elimination-partial least square (UVE-PLS) were used to develop the statistical models and to generate the steric and electrostatic contour plots (Edache et al. 2022). The resulting models were then cross-validated using the leave one out (LOO), leave two out (LTO), and leave many out (LMO).

Molecular docking screening
A molecular docking simulation was conducted between the thirty-six arylimidamide-azole analogs, the reference drug (pentamidine), and the pyridoxal kinase (PdxK) receptor (PDB: 6K91) using Molegro Virtual Docker (MVD). PdxK in 3-D form was retrieved from the protein data bank and prepared on the Molegro Virtual Docker by removing water molecules, ligands, and cofactors associated with the protein structure. The software allows for the repair (rebuild) of all affected residues. The receptor's binding cavities were defined and that which has the largest volume and surface (volume: 398.336 and surface: 952.32) was adopted for the docking. All ligands were imported in PDB file format and prepared appropriately. The simulation was carried out using the parameter settings available in Table 3. The predicted ligand-protein interaction profiles were then visualized using the Biovia Discovery Studio Visualizer. A similar method was earlier reported elsewhere (Ibrahim et al. 2021;Abdullahi et al. 2022).

Evaluation of pharmacokinetic properties
The prediction of pharmacokinetic properties plays a very critical role in the early stage of drug discovery. This is because only molecules which demonstrate good ADMET and drug-likeness properties reach the preclinical research phase (Ugbe et al. 2022a, b). Therefore, six arylimidamide-azole hybrids (21, 22, 26, 31, 33, and 36) with the highest anti-leishmanial activities were subjected to drug-likeness and ADMET tests using two online web servers; http:// www. swiss adme. ch/ index. php and http:// biosig. unime lb. edu. au/ pkcsm, respectively. Similar tests were performed on the newly designed compounds (36a-36 g) to ascertain their drug-likeness properties. Lipinski's RO5 also called the Pfizer rule is a well-established provision determining the oral bioavailability of a given compound (Lipinski et al. 2001;Lawal et al. 2021). Hence, the various compounds would be subjected to the RO5 criterion to ascertain their oral bioavailability.

Maximum iteration 1500
Maximum population size 50

Pose generation
Energy threshold 100 Simplex evolution

Maximum steps 300
Neighbor distance factor 1.00

Multiple poses
The maximum number of poses returned 5 Enable energy threshold 0.00 Cluster similar poses RMSD threshold: 1.00

Ligand-based drug design (LBDD)
QSAR has over the years been utilized for ascertaining the relationships between the bioactivities of chemical entities and their physicochemical properties, to obtain a reliable model which could be used for predicting the antiproliferative activities of new chemical compounds. This is governed by the variation of structural properties across compounds which in turn is responsible for the difference in their biological activities. 3-D QSAR has emerged as a more sophisticated tool for designing new molecules than the classical QSAR study (2-D) because it exploits the 3-D properties of the ligands to predict their activities (Verma et al. 2010). In this study, seven compounds (36a-36 g) were designed via LBDD using the template compound (36) majorly by structural adjustments involving the addition and/or removal of certain chemical moieties or substituent groups based on the information provided by the steric and electrostatic fields' contour maps. The molecular structures of the new molecules were prepared according to the procedure earlier described under the 'molecular geometry minimization' section. Molecular docking was conducted between the newly designed compounds and the target receptor (PdxK), using MVD, while utilizing the Biovia Discovery Studio to visualize the resulting pharmacological interactions. The bioactivities of the new compounds were predicted by the built 2-D QSAR model equation.

2-D QSAR modeling
Theoretical modeling of thirty-six arylimidamide-azole derivatives was conducted to establish a quantitative relationship between their structures and their inhibitory activities. As a result, a five-descriptor QSAR model was built (Eq. 13) with the descriptors well described in Table 4.  The outcome of the internal and external validation assessment conducted on the built model is available in Table 5.
The computed descriptors, observed activities (pIC 50 ), and the predicted activities together with their residuals are presented in Table 6. Also, a plot of predicted activities versus experimental activities for the training set and test set is shown in Fig. 1, while Fig. 2 shows the plot of standardized residuals against experimental activities (pIC 50 ).
Furthermore, the result of Pearson's correlation statistical analyses performed on all five descriptors' values in the model building set is reported in Table 7. Another powerful validation method called Y-randomization or Y-scrambling (13) test was equally applied, and the result is presented in Table 8. Finally, Fig. 3 shows the plot of standardized residuals against the leverages otherwise known as William's plot determined to show the model's applicability domain space.

3-D QSAR modeling (CoMFA)
Molecular structural alignment represents a crucial factor in determining the predictive capacity of a built 3-D QSAR model. Figure 4a, b shows the molecular structure of the alignment template (compound 32) and the aligned structures as obtained from the atom-based superimposition of the remaining structures on the template. Two models were developed: each from the FFD-SEL and UVEPLS approaches. Some significant statistical parameters calculated for both the 3-D models are presented in Table 9. Presented in Table 10 are the O3A scores, experimental pIC 50 , and the predicted pIC 50 available from the two models, as well as their residual values. Additionally, a plot showing the correlation between predicted and experimental activities for training and test sets was obtained for both models and is presented in Fig. 5a, b. The CoMFA model equation (FFDSEL) is summarized graphically as 3-D contour maps as shown in Figs. 6a, b and 7a, b.

Molecular docking study
The results (MolDock scores) of the molecular docking simulation conducted between the receptor (PdxK) and the 36 arylimidamide-azole hybrids as well as the reference drug (pentamidine) using the Molegro Virtual Docker are earlier reported in Table 1. The predicted pharmacological interactions of compound 36 (template) and pentamidine, with the target receptor, are shown in Fig. 8.

Pharmacokinetics study
The results of the drug-likeness and ADMET properties investigation conducted on the selected  Tables 11 and 12.

Ligand-based drug design (LBDD)
A total of seven compounds were designed using the ligand-based drug design approach. A combined 2-D and 3-D QSAR approach was adopted, where information provided by the 3-D QSAR (CoMFA) contour maps was used for the design, while the inhibitory activities of these new compounds were predicted using the 2D-QSAR model equation. Their molecular structures predicted activities and docking scores are presented in Table 13, while Table 14 shows their predicted pharmacological interaction profiles with PdxK. Also, the binding interactions of these compounds (excluding 36f and 36 g) with the active sites of PdxK as visualized via Biovia discovery studio are presented in Figs. 9, 10, 11. Their predicted drug-likeness and ADMET properties are earlier reported in Tables 11 and 12.

2-D QSAR modeling
The values of VIF reported in Table 7 range from 1 to 5 for all 5 descriptors, indicating that the built QSAR model was statistically substantial and is therefore stable and acceptable. Each descriptor has an absolute t-statistics value of greater than 2, which shows that the selected descriptors were significant (Adeniji et al. 2018). Additionally, the calculated p values (Table 7) of these descriptors in the model at a 95% confidence level were less than 0.05. This is therefore in support of an alternative hypothesis that states that there is a relationship between the descriptors in the model equation and the compounds' inhibitory activities at p ˂ 0.05. The strength and direction of each descriptor's contributions in the built QSAR model are determined by the value of the mean effect (ME) ( Table 4). The magnitude of ME signifies the extent to which the descriptor influences the molecule's inhibitory activity, while the sign of ME shows the direction of influence. As shown in Table 4, SpMin8_Bhv and AATS8p have the largest positive ME, which implies elevating their values will mean a corresponding increase in the molecule's inhibitory activities. ATSC3c, VR1_Dzv, and MDEC-22 on the other hand have negative ME values, signifying that an increase in these descriptors' values will lead to a decrease in the molecule's antiproliferative activities.
The application of GFA in conjunction with the MLR approach led to the development of the 2-D QSAR model with a total of five descriptors as shown in Eq. 13. The built model was adjudged to satisfy the requirement of a good QSAR model. The low residual values between the experimental and predicted activity values for both training and test set compounds as shown in Table 6 indicate a high predictive strength of the built model. The R 2 values of 0.9614 for the training set and 0.6766 for the test set from Fig. 1 were similar to those obtained from GFA (0.9615 and 0.6766) and MLRplusValidation analysis (0.9615 and 0.6766) available in Table 5, which means that the R 2 values were computed accurately. Figure 1 shows how well the predicted activities correlated with the experimental activities as suggested by the grouping of points along the line of best fits, implying that the built model is very reliable and robust in predicting the activities of new chemical entities. A random distribution of standard residuals on both sides of zero as observed in Fig. 2 shows that the built model is not associated with any systematic error.
Furthermore, the result of the Pearson correlation analysis available in Table 7 shows that no pair of descriptors is significantly inter-correlated as suggested by the low values of the correlation coefficients. This was also observed by Adeniji et al. (2018).
The values of R 2 and Q 2 obtained from the y-randomization test (Table 8) were significantly low, and this implies that the developed model is stable, robust, and reliable. The Y-randomization coefficient, c R 2 p (0.877748) greater than 0.50, signifies that the built model is powerful and not inferred by chance. The plot of standardized residuals against leverages (William's plot)    Fig. 3) shows that all the compounds were contained within the ± 3.0 square area of standardized cross-validated residual, which clearly shows that no outlier is present in the data set. However, compounds 1 and 32 have leverages greater than the computed warning leverage (h* = 0.72) and were said to be influential molecules.

3-D QSAR modeling
The 3-D QSAR modeling began with the atom-based alignment of structures using compound 32 as the alignment template. The alignment process involves an early step that provided all the 36 compounds the opportunity of being selected as the alignment template based on the compound with the highest Open3DAlign (O3A) score. The O3A scores of the various compounds are included in Table 10. Compound 32 was found with the highest O3A score of 6594.08 and hence selected as the template upon which the remaining structures were superimposed. Two CoMFA models were built using the FFDSEL and UVEPLS approaches, with their computed statistical parameters tabulated in Table 9. The parameters were computed for five principal components (PCs) among which the fifth PC (PC = 5) performed relatively better. The statistical parameters available in Table 9 were those associated with PC 5. The predictive strength of the regression models on new datasets of compounds can be estimated by cross-validation (Grohmann and Schindler 2008). A cross-validated coefficient of correlation (Q 2 ) ≥ 0.50 indicates a good QSAR model. Here, three types of Q 2 were calculated; leave one out (LOO), leave two out (LTO), and leave many out (LMO), together with their associated standard error of prediction (SDEP). Although both models passed the cross-validation test (i.e., Q 2 > 0.5), relatively higher values of correlation coefficient (R 2 ), Q 2 LOO, Q 2 LTO, Q 2 LMO, and the Fischer's statistics (F-test) were associated with the UVEPLS CoMFA model. Also, the calculated standard error of correlation (SDEC) and SDEP were lower for UVEPLS than for FFDSEL. A linear correlation between the CoMFA descriptors (independent variables) and the activity values (dependent variables) was established by the PLS analysis method. The lower residual values between the predicted and observed activity values (Table 10) show a strong predictive strength of both models. This was supported by the clustering of points along the lines of best fits in the plots of predicted pIC 50 against the experimental pIC 50 for both models (Fig. 5a, b).
The CoMFA QSAR equation is summarized graphically as a 3-D contour map, which shows the regions within the molecules' 3-D structural space where steric and electrostatic fields are associated with extreme values. The underlying principle behind CoMFA is that variations in the shape and strength of non-covalent interaction fields surrounding the molecules, such as steric or electrostatic fields, can be related to changes in binding affinities (Kakarla et al. 2016). Therefore, molecular fields are key factors in binding affinity.
The contour maps presented in Figs. 6, 7 were those generated from the UVEPLS CoMFA model for compound 36. The steric and electrostatic field contributions are 68.19% and 31.81% for FFDSEL, and 68.96% and 31.04% for UVEPLS, respectively (Table 9), which has only a very insignificant difference between both models.   From the steric field contour maps available in Fig. 6a, b, the red contours represent the regions of unfavorable steric bulk, while blue contours represent the regions of favorable steric bulk. An identified region of unfavorable steric bulk in compound 36 is the number 2-position of the alkyl chain linking the biphenyl group and the imidazole group. This means that if the hydrogen in that position is replaced with a more bulky group like methyl, ethyl, isopropyl, etc., the compound's activity or binding affinity will decrease. On the other hand, more steric bulk favorable regions were identified, which include the alkoxy substituents on the biphenyl ring system and the 4-position of the pyridine ring system. This implies that increasing the steric bulk around those areas will enhance the inhibitory activity or binding affinity of the molecule. From the electrostatic field contour maps available in Fig. 7a, b, yellow contours represent regions favored by high electron density or unfavorable to electron-withdrawing substituents, while green contours represent regions of unfavorable high electron density or favorable to electron-withdrawing groups. Two regions of low electron density were identified to include the number 4-position of the alkyl chain linking the biphenyl group to the imidazole group, and the number 2-position of the phenyl ring which links to the amide group. These regions need not be too electron-dense, and hence, the introduction of electron-withdrawing groups at these positions will enhance the molecule's inhibitory activity or binding affinity. Regions of high electron density were identified around positions 2 and 3 of the phenyl ring which links to the alkyl chain via the oxygen, and positions 3 and 5 of the pyridine ring system. These regions need not be less electron-dense, and hence, electron-donating groups will keep these regions at a high electron density which in turn will enhance the molecule's inhibitory activity or binding affinity.

Molecular docking study
The binding affinity of protein-ligand interaction is necessary to describe how well drugs bind to the target receptor. The spontaneity of the binding process and how well drugs can fit into the receptor's binding cavities to form the most stable receptor is shown by the negative value of the binding affinity (Ugbe et al. 2022a, b). Here, the results (binding affinities) obtained from the docking study using Molegro Virtual Docker are expressed as MolDock scores in kcal/mol as included in Table 1. On average, the result shows that binding affinity correlates well with inhibitory activity. The screening identified compounds 1 and 36 as showing the least and highest binding affinities of -142.632 kcal/mol and -198.260 kcal/mol, respectively. Also, the predicted binding affinity (MolDock score = − 141.793 kcal/mol) of the reference drug (pentamidine) was lower than those of the arylimidamide-azole hybrids, an indication that the various compounds bind better with the active sites of the protein (PdxK) than pentamidine. This claim can further be supported by comparing the binding interaction profile of 36 with that of pentamidine as presented in Fig. 8 and Table 14. More interactions were visible from the binding profile of compound 36 with PdxK. For 36_PdxK, three conventional hydrogen bonds were formed involving GLN-258 (interaction distance of 2.78 Å), THR-229 (1.95 Å), and SER-47 (2.42 Å), Two carbon-hydrogen (C-H) bonds were visible with HIS-46 and ASP-124, a π-donor hydrogen bond with GLY-230, π-anion electrostatic interaction with ASP-124, and several hydrophobic interactions. Pentamidine_PdxK on the other hand shows two conventional hydrogen bonds with THR-229 (2.53 Å) and SER-47 (2.52 Å), one C-H bond involving TYR-129, a π-donor hydrogen bond with GLY-228, and one π-alkyl hydrophobic interaction with VAL-19.

Pharmacokinetics properties prediction
Lipinski's RO5 for oral bioavailability has a broad application in the discovery of new drug molecules (Ugbe et al. 2022a, b). It asserts that a drug molecule may likely not be orally bioavailable when it has hydrogen bond donors (HBDs) of greater than 5, hydrogen bond acceptors (HBAs) > 10, molecular weight (MW) > 500, and lipophilicity (MLOGP > 4.15 or WLOGP > 5) (Lipinski et al. 2001). Whenever a molecule passed at least three of the four provisions of the RO5, it is said to obey Lipinski's rule for oral bioavailability (Lawal et al. 2021). From Table 11, it was observed that all the selected arylimidamide-azole hybrids passed the drug-likeness test (Lipinski RO5). Five of the seven newly designed analogs pass the Lipinski criteria. Only 36f and 36 g failed the drug-likeness test  Table 12 Predicted ADMET properties of some selected arylimidamide-azole derivatives and the newly designed analogues BBB-Blood-brain barrier, CNS-central nervous system, HIA-human intestinal absorption, Skin-skin permeability, LogBB-the logarithmic ratio of brain to plasma drug concentration, LogPS-blood-brain permeability-surface area product, CYP-34A-cytochrome p450 isoform, CYP-2D6-cytochrome p450 isoform, Ssubstrate, I-inhibitor, MRTD-maximum recommended tolerated dose, TCE-total clearance  The predicted ADMET properties available in Table 12 show a very high human intestinal absorption for all tested compounds with 36 and the new compounds reaching 100%. All the tested molecules are substrates and inhibitors of p-glycoprotein, indicating that these molecules may easily mediate to reach their target sites. Skin permeability is a key factor in transdermal drug delivery development. Values of skin permeation constant LogKp > -2.50 connotes low skin permeability. As a result, the various compounds tested showed LogKp values < -2.50, indicating good skin permeability. Drug molecule penetration through the blood-brain barrier (BBB) and central nervous system (CNS) comes with certain criteria. For a drug molecule to be able to pass through the BBB and CNS readily, the logarithmic ratio of brain to plasma drug concentration (logBB) must be > 0.3 and the blood-brain permeability-surface area product (logPS) be > -2, respectively. As a result, all the selected arylimidamide-azole hybrids do not readily penetrate the BBB and the CNS. On the other hand, all the newly designed compounds are poorly distributed to the brain but are said to readily penetrate the CNS except 36b. Furthermore, some group of enzymes called cytochrome P450 enzymes are important in the body to facilitate drug metabolism and to help in their excretion. The two major isoforms enhancing drug   Table 12 is the maximum recommended tolerated dose (MRTD) predicted for the various molecules. MRTD value of ≤ 0.477 log (mg/kg/day) is considered low, while a value > 0.477 log (mg/kg/day) is considered high. The overall druglikeness and ADMET properties of the design template (36) and the newly designed compounds (except for 36f and 36 g) showed good pharmacokinetic profiles. Therefore, these molecules could be considered potential drug candidates for the treatment of leishmaniasis.

Ligand-based drug design (LBDD)
As stated earlier, the information encoded in the 3-D QSAR contour maps of the design template (36) was utilized in designing seven new arylimidamide-azole derivatives (36a-36 g). 36a was obtained by replacing the isopropoxy groups of the template with a bulky tert-butoxy group to increase the steric bulk as observed from the contour maps. 36b was designed by reducing the steric bulk around the 2-position of the alkyl (hexyl) chain with the introduction of two double bonds (conjugated) to take off four hydrogens around that position. 36c reduces the two double bonds in 36b to one. 36d replaces the isopropoxy groups in 36c with tert-butoxy groups. 36e introduces the methyl group (a more bulky group than hydrogen) at the para-position of the pyridine ring system and also takes into account the electrostatic field effect by introducing the chloro group (electron-withdrawing group) at that position. 36f further introduced the chloro group at the 4-position of the hexyl chain to deplete the electron density in that region. Finally, 36 g further introduced two methyl groups at positions 3 and 5of the pyridine ring system. The results of activity prediction and molecular docking conducted on the newly designed compounds (Table 13) showed that all the compounds have higher predicted inhibitory activities than the template (36) with only 36b and 36 g showing lower MolDock scores than 36. Compound 36e showed the highest MolDock score of -208. 595 kcal/mol and the second-highest predicted pIC 50 of 7.375 only behind 36f with a predicted pIC 50 value of 7.670. Also, 36f and 36 g failed the drug-likeness test as both violated 2 out of 4 provisions of the Lipinski RO5. Hence, the interaction profile of compound 36e will be discussed further. One C-H bond involving ASP-125 at a distance of 2.69Å and two π-donor hydrogen bonds with GLY-228 (3.32Å) and GLY-230 (2.83Å) were visible in the binding interaction profile of 36e with PdxK. Electrostatic interaction types are attractive charge with ASP-231 and π-anion with HIS-46 and ASP-124. Others are hydrophobic interactions such as π-π stacked with HIS-46, alkyl with VAL-121, π-alkyl with HIS-46 and TYR-12, and halogen interaction with ASP-231. Interestingly, unlike the template and reference drug (pentamidine), most of the newly designed compounds have no unfavorable clashes with some amino acid residues, indicating a possible elimination or lowering of side effects that could be associated with the newly designed hybrids. Therefore, compound 36e and the other newly designed analogs (except 36f and 36 g) have shown good pharmacological and pharmacokinetic properties and are hence worthy of selection for further practical evaluation in the laboratory as pyridoxal kinase inhibitors.

Conclusions
In this study, a five-descriptor 2-D QSAR model and a 3-D QSAR (CoMFA) model were developed with 36 arylimidamide-azole hybrids, both of which were found to satisfy the requirement for internal and external validation tests. The anti-leishmanial activities of the various compounds were well predicted by both models. A combined 2-D and 3-D QSAR approach was utilized to design and predict the inhibitory activities of seven new arylimidamide-azole analogs using compound 36 as the template. The molecular docking screening conducted between the 36 arylimidamide-azole compounds and the target receptor, pyridoxal kinase (PdxK), revealed compounds 1 and 36 as having the least and highest Mol-Dock scores of − 142.632 kcal/mol and − 198.260 kcal/ mol, respectively, both higher than MolDock score of -141.793 kcal/mol reported for the reference drug pentamidine. The new compounds bind excellently into PdxK's cavities with binding affinities (MolDock scores) in the order: 36e (− 208.595 kcal/mol) > 36f > 36a > 3 6c > 36d > 36 (Template) > 36b > 36 g (− 187.155 kcal/ mol), while the predicted pIC 50 follows the order: 36f (7.670) > 36e > 36 g > 36c > 36d > 36a > 36b > Template (36) (6.344). Also, the newly designed analogs showed a good pharmacokinetic profile and obeyed Lipinski's RO5 for oral bioavailability except for 36f and 36 g. Special emphasis on 36e because it appears to be the most consistent with the various employed validation protocols, being that it possessed the highest binding score, second-highest predicted pIC 50 , showed excellent pharmacokinetic properties and binds adequately with the target protein. Therefore, this work has provided useful information on the arylimidamide-azole hybrids as antileishmanial agents.