Drug design of new anti-EBOV inhibitors: QSAR, homology modeling, molecular docking and molecular dynamics studies

2.1 Data collection and preparation

A dataset of 86 amodiaquine derivatives with EBOV-GP inhibitory activity were collected from the published work by Sakurai et al. (Sakurai et al., 2018). After drawing the chemical structures of these derivatives using ChemDraw software (PerkinElmerInformatics, Inc, ChemDraw, 2021), their geometries were optimized in the gas phase using the AM1 semi-empirical method. The optimization procedure was performed using the Gaussian 09 software package (Gaussian 09, n.d.). The optimized molecular structures were converted into isomeric SMILES notations for the modeling process of the CORAL model using the Openbabel program (O’Boyle et al., 2011). For the GA-MLR model, the molecular descriptors of the optimized 3D structures were calculated using AlvaDesc (Mauri and Bertola, 2022). To avoid multicollinearity in the QSAR model, we reduced the number of generated variables by removing descriptors with a correlation coefficient greater than 0.95 and retaining only one of the pairs. As a result, a total of 832 descriptors were selected from the original set of 4149 descriptors. The experimental data value of each molecule (half-maximal inhibitory concentration: IC₅₀) was converted into negative logarithmic scale (pIC₅₀). These converted data were then utilized as variables for building QSAR models. The molecular structures and corresponding IC₅₀ data are listed in Table S1 and their SMILES notation and converted pIC₅₀ can be found in Table S2.

2.2 CORAL QSAR model

CORAL QSAR modeling software (Toropova and Toropov, 2014) utilizes SMILES attributes for endpoint prediction, based on optimal descriptors, namely correlation weights (CWs), along with the balance-of-correlation approach. Four splits were randomly made from the entire dataset with each split in turn divided into four partitions: Training (Tr: 35 %), Invisible Training (Inv-Tr: 35 %), Calibration (Calib: 15 %), and Validation (Valid: 15 %). The specific functions of each set in the development of QSAR models are described in detail in the literature (Toropova and Toropov, 2014; Oubahmane et al., 2023). Table 1 illustrates the distribution and percentage of the identity of the compounds in the four splits.

The development of QSAR models based on SMILES notation uses the equation of the optimal descriptor of correlation weights (Eq.1).

Following the calculation of all CWs, the linear regression technique was then employed to develop QSAR models, as indicated in equation (Eq.2).

In order to optimize the QSAR modeling, based on the Monte Carlo method, two different target functions, namely TF₁ (without the IIC) and TF₂ (with value W_IIC = 0.3) were used to generate the models. TF₁ was calculated using the equilibrium of the correlation method (according to equation (Eq.3)), and the Index of Ideality of Correlation (IIC) (Liman et al., 2022; Toropov et al., 2019; Toropova and Toropov, 2017) was combined with TF₁ to formulate the modified target function TF₂ (Eq.4).

The IIC, given in equation (Eq.5), was proposed as a measure for evaluating the predictive ability of the developed QSAR models. In particular, it improves the precision of the model, which is determined by the determination coefficient (R²) and the mean absolute error (MAE). The W_IIC coefficient, which represents the IIC weight, can change the extent of the influence of the IIC on the Monte Carlo optimization process. The optimal value of W_IIC is influenced by two aspects: the diversity of the molecules and the characteristics of the endpoint.

Where;

Δ_k measures the quality prediction of the k^th molecule.

To construct robust QSAR models, the optimal threshold (T*) and the number of epochs (N*) had to be determined by evaluating the most favorable statistical metrics for the calibration set. In the search for the ideal T* and N*, ranges from 1 to 10 were used for the threshold and 1 to 30 for the number of epochs in the optimization process, using three optimization probes. For this study, the W_IIC was set to 0 for TF₁ and adjusted to 0.3 for TF₂.

2.3 GA-MLR QSAR model

The selection of the most relevant descriptors, from the entire computed ones, is a crucial step in QSAR modeling. This selection was achieved by stepwise MLR technique via XLSTAT (Addinsoft, 2020). A total of 832 molecular descriptors were filtered out to get a final selection of five relevant molecular descriptors. Using these descriptors, the MLR method with the ordinary least squares (OLS) algorithm, implemented in the QSARINS software (Addinsoft, 2020; Gramatica et al., 2013) was applied to establish a linear correlation between the pIC₅₀ endpoints of the molecules and their corresponding molecular descriptors. The dataset was divided into a training set with 60 molecules and a test set consisting of 26 molecules, with a distribution ratio of 70 % and 30 % accordingly. The GA-MLR models were developed using the following parameters: all subsets up to a maximum of 5, a maximum generation of 10 000, and a mutation probability of 0.05. The remaining settings were configured with their default values.

2.4 QSAR model validation

The validation process is a crucial stage in QSAR establishment to assess the accuracy of the model and to make reliable predictions about the activity of new molecules. This process is considered the most important aspect in testing the robustness, predictability, and reliability of a QSAR model. Four steps are typically used to validate the generated model, including (a) internal validation or cross-validation using the training set, (b) Y-randomization, (c) independent validation using the test set, and (d) applicability domain (AD) evaluation. The validation procedures and criteria for the CORAL-based QSAR and GA-MLR QSAR models are well discussed in our earlier works (Gramatica et al., 2014; Hdoufane et al., 2019, 2022; Oubahmane et al., 2023; Toropova and Toropov, 2014).

2.5 EBOV-GP structure preparation

The experimental 3D structure of EBOV-GP (PDB ID: 5JQ7) has some missing residues that are essential for its functional mechanism in three regions (LYS191-SER210, THR284-GLY286 and LYS294-ARG302). To fill these missing residues and facilitate subsequent molecular docking and MD simulations, we used the SWISS-MODEL web server for homology modeling to reconstruct the full structure of EBOV-GP (Schwede et al., 2003). The server applies a rigorous protocol that includes the use of BLAST and HHblits to find structural templates from its extensive library, the SWISS-MODEL Template Library (SMTL), which is sourced from the PDB. This ensures that high quality templates are provided, as determined by the Global Model Quality Estimate (GMQE) and the Quaternary Structure Quality Estimate (QSQE). After selecting the most suitable template, which in this case was the same PDB ID due to its high similarity to the target structure, the incomplete EBOV-GP structure was aligned with the template. Subsequently, the missing loop regions were carefully filled in using OpenStructure and the ProMod3 modeling algorithm (Studer et al., 2021). The accuracy and reliability of the final EBOV-GP model was checked using the Qualitative Model Energy ANalysis (QMEAN) assessment function to ensure that the final model was of high quality and suitable for the intended advanced computational analyses (Benkert et al., 2011).

2.6 Molecular docking

Molecular docking analysis was performed to investigate the ligand–protein interactions of the designed compounds and the reference lead compound (78) within the binding site of the modeled EBOV-GP protein. Prior to the docking process, the co-crystallized ligand (Toremifene: TOR) was used to generate the grid-box coordinates of the binding site of the EBOV-GP protein. The grid box coordinates are as follows: xyz coordinates of −45, 16.75, −9, and grid box dimensions of 24 ų. The protein was then prepared by removing water molecules and the co-crystallized ligand and adding polar hydrogens and Gasteiger charges using AutodockTools (Morris et al., 2009). The protein and compounds were saved in PDBQT format using the same software. Molecular docking analysis were performed using AutoDock Vina software (Trott and Olson, 2010) with exhaustiveness of 20. The docking results were analyzed by the Discovery Studio visualizer and Pymol softwares (BIOVIA Discovery Studio Visualizer, n.d.; Yuan et al., 2017).

2.7 Molecular dynamics simulations

The apo form of the EBOV-GP protein and the designed compounds complexed with the EBOV-GP protein were subjected to all-atom molecular dynamic simulations using GROMACS 2021.3 software (Van Der Spoel et al., 2005). Prior to the MD simulations, the CHARMM-GUI web server (Jo et al., 2008) was utilized to generate initial input parameters, implementing the CHARMM36 force field (Huang et al., 2017). Each complex was solvated in TIP3P water within a rectangular grid box (Jorgensen et al., 1983) and counterions (Na⁺, Cl^-) were added to maintain a 0.15M salt concentration through Monte Carlo ion displacement. The size and shape of the systems were considered when the periodic boundary conditions (PBC) were imposed. Non-bound interactions were addressed using Lennard-Jones and Coulomb potentials with a cut-off distance of 12 Å, and the Verlet cut-off strategy was used to buffer the neighbor search list. The Particle-Mesh Ewald (PME) method (Darden et al., 1993) was used to address the long-range electrostatic interactions. The novel Linear Constraint Solver (LINCS) algorithm (Hess et al., 1997) constrained all covalent bonds, including hydrogen atoms. Energy minimization of each system was carried out using the steepest descent algorithm with a maximum of 50,000 steps and a maximum force of 10.0 KJ/mol. The systems were equilibrated using a two-step process: first, a canonical (NVT) ensemble was employed to equilibrate the temperature at 303 K for 500 ps, followed by an isothermal-isobaric (NPT) ensemble for another 500 ps to equilibrate the pressure, both set at 303 K and 1.01325 bar respectively. The Nose-Hoover thermostat (Hess et al., 1997; Nosé, 1984) and Parrinello-Rahman barostat (Parrinello and Rahman, 1981) were used to regulate the temperature and pressure respectively.

The production run was conducted for 100 ns to analyze the structural stability of the designed molecules within the binding site of EBOV-GP protein. The structural stability of the designed molecules within the EBOV-GP was assessed using the root mean square deviation (RMSD), the radius of gyration (Rg), the root mean square flexibility (RMSF), and the solvent accessibility area (SASA) based on the dynamics trajectory results. The data obtained were visualized using Xmgrace software (GRACE - GRaphing, Advanced Computation and Exploration of data, (n.d.). https://plasma-gate.weizmann.ac.il/Grace/ (accessed March 4, 2024).

2.8 Drug likeness and ADMET prediction

The evaluation of ADMET has a great importance in the initial phases of pharmaceutical exploration. It is expected that a therapeutic drug of better quality would demonstrate exceptional activity at the target receptor, as well as better ADMET profile at therapeutic dosage. Therefore, it is fundamental to assess the pharmacokinetic properties of hit compounds to minimize the risk of subsequent medication failure. The qualities related to drug-likeness are important to understand the distribution of a molecule in an organism, which in turn affects its pharmacological effectiveness. ADMET predictions for the designed compounds were performed via the AdmetLab2, ProTox-II, and Osiris Property Explorer servers (Xiong et al., 2021; Banerjee et al., 2018; Property explorer, 2024).

3.1 CORAL-based QSAR model

A total of eight QSAR models were created by performing four random splits. The two different target functions, namely TF₁ and TF₂, were used to generate these models. The obtained statistical parameters from the analysis of these QSAR models suggest that the inclusion of W_IIC = 0.3 significantly increases the influence of IIC in the Monte Carlo optimization process. Table 2 presents the statistical parameters for all splits, indicating the robustness and validity of the models according to the criteria established by Tropsha et al. (Golbraikh and Tropsha, 2002) and Ojha et al. (Ojha et al., 2011). Fig. 1 displays the correlation between the experimental and predicted pIC₅₀ values across the four splits. The QSAR model derived from split 3 was chosen as the optimal model based on its superior validation metrics ( $R^2$ = 0.9129; $Q^2=0.8848$ ). The mathematical model equation for split 3 is given below (Eq.9):

Close

Fig. 1

Experimental vs. calculated pIC₅₀ values for the four Splits.

Export to PPT

3.2 GA-MLR QSAR model

The GA-MLR model was built based on five selected descriptors from a comprehensive set of descriptors, namely IC2, SpMin7_Bh(s), P_VSA_Charge_7, CATS2D_05_DD, and T(O..O). The description of each descriptor was given in Table S5. These descriptors were found to significantly explain the structure activity relationship (SAR) between the pIC₅₀ activity and the studied compounds and were subsequently used to construct the QSAR model. The generated model using the GA-MLR approach (Eq.10) and its corresponding statistical parameters are shown below:

• N_Tr = 60, R² _Tr = 0.8490, Q²_LOO = 0.8148

• N_ext = 26, R²_ext = 0.8661, MAE_ext = 0.2135.

• Q²_F1 = 0.8681, Q²_F2 = 0.8599, Q²_F3 = 0.8840, CCC_ext = 0.9229, F = 77.6036, s = 0.3138.

N_Tr and N_ext indicate the total number of samples used for the training and the test sets, respectively. CCCext is the external concordance correlation coefficient, while F refers to F-statistic. Additionally, the standard deviation is represented by “s” and the external validation criteria are expressed through $Q_{F1}^2$ , $Q_{F2}^2$ and $Q_{F3}^2$ (Chirico and Gramatica, 2012; Hess et al., 1997; Nosé, 1984).

The evaluation of the indicated parameters in the established GA-MLR model satisfies the validation requirements outlined in the OECD guidelines. Moreover, the scatter plot shown in Fig. 2a displays a good correlation between the experimental and predicted pIC₅₀ endpoints. To confirm the reliability of the built model, the AD was evaluated using the leverage method, and presented through the Williams plot in Fig. 2b. The graphic includes dashed lines to indicate the cutoff value of ±3 s.d., and the warning line for the outliers (h*) was found to be 0.29. The findings from the Williams plot provide definitive evidence that the AD covers all compounds, except for compound 65. This observation can be attributed to the low activity of compound 65, as evidenced by its high value (IC₅₀ = 20 μM), which is 69-fold higher than that of the lead compound (IC₅₀ = 0.29 μM).

Close

Fig. 2

(a) Experimental vs. predicted pIC₅₀ values; (b) William’s plot.

Export to PPT

3.3 Mechanistic interpretation

The integration of mechanistic interpretations is a pivotal component within the framework of the OECD. These interpretations are helpful for identifying and evaluating the chemical features that contribute to the increase or decrease of a defined endpoint value. To provide a comprehensive understanding of the mechanistic interpretation of the CORAL model, numerous iterations of Monte Carlo optimization were used. The chemical features, obtained from SMILES descriptors, display positive (CWs) and recognize factors that support an increase in the pIC₅₀ value. In contrast, negative CWs identify factors that encourage a reduction in the pIC₅₀ value. While SMILES qualities that include both positive and negative CWs values remain undefined. Table S6 presents a complete summary of the main factors that contribute to the increase or decrease of pIC₅₀ values, associated with their corresponding correlation weights (CWs). These results are derived from the three separate runs of the established QSAR model for split 3.

Based on the SAR study of Sakurai et al. (Sakurai et al., 2018) and on the insights obtained from Table S6, the factors increasing the inhibitory activity were incorporated into the lead compound 78 along with compounds 80 and 83. These two molecules exhibit low cytotoxicity and high pIC₅₀ values and are closely aligned with the lead compound. While the promoters associated with a decrease were omitted. The influential promoters were systematically investigated at three distinct sites (R₁, R₂, R₃) within the lead compound to strategically design novel inhibitors targeting the EBOV-GP, as depicted in Fig. 3. Therefore, 26 novel molecules were designed (Table S7). Table S8 gives the calculated activities of the newly designed molecules by the two QSAR models. The CORAL model predicted a pIC₅₀ range [7.1239–7.9595] for all the designed compounds. The GA-MLR QSAR model confirmed the inhibitory potential of the newly developed compounds.

Close

Fig. 3

Scaffold chemical structure derived from the lead compound with highlighted substitution sites (R₁, R₂, R₃) maintaining crucial Iodine for inhibition activity.

Export to PPT

3.4 Modeling missing residues as loop regions

The loops are involved in many biological functions and constitute key components of protein structure. The loop modeling of the EBOV-GP (PDB ID: 5JQ7), successfully filled the missing residues in three critical regions (LYS191-SER210, THR284-GLY286 and LYS294-ARG302). These regions are essential for the functional mechanism of the EBOV-GP as depicted in Fig. 4. The process involved aligning the incomplete structure of EBOV-GP with its corresponding higher-resolution template, and the missing residues were completed using the ProMod3 modeling algorithm. The accuracy and reliability of the reconstructed EBOV-GP was confirmed by QMEAN and GMQE values of 0.69 and 0.85, respectively. Importantly, the loop model of the missing residues in region 1 is positioned and oriented upwards to allow the ligands to access the binding site. This structural arrangement is in agreement with the findings of previous research (Ren et al., 2018), suggesting that the modeled loop does not obstruct the binding pocket and thus facilitates potential inhibitor interactions. These findings indicate that the completed EBOV-GP structure is well constructed, confirming its suitability for subsequent molecular docking and MD simulations.

Close

Fig. 4

(a) The alignment of the uncompleted structure of EBOV-GP (PDB ID: 5JQ7, light pink) and the modeled EBOV-GP (light blue) containing the co-crystalized ligand (TOR). The modeled structure of the missing residues in the three regions: (b) Region 1 (LYS191-SER210, blue), (c) Region 2 (THR284-GLY286, magenta) and Region 3 (LYS294-ARG302, orange).

Export to PPT

3.5 Molecular docking

To validate the efficacy of our design approach, we employed molecular docking methodology to explore the ligand–protein interactions involving the newly designed compounds, the reference lead compound, and the co-crystallized ligand (TOR) with the modelled EBOV-GP protein. A crucial step in this validation process was the re-docking analysis of the TOR ligand to assess the accuracy of the docking protocol. The RMSD between the x-ray and predicted poses of TOR was determined to be 0.664 Å (Figure S1), confirming the appropriateness of the docking protocol in reproducing native poses within the acceptable range of less than 2 Å. Subsequently, a comprehensive docking study was conducted for all designed compounds against the modelled EBOV-GP, and the resulting binding affinity values are detailed in Table S8. Notably, ten designed molecules exhibited higher binding affinity values compared to both the co-crystallized ligand (TOR: −8.4 Kcal/mol) and the lead compound (78: −7.4 Kcal/mol). The chemical structures of the top five newly designed compounds (D1-D5) ranked based on their binding affinities are presented in Fig. 5. While, Fig. 6 illustrates their 2D and 3D representations within the binding site of EBOV-GP. Their binding affinity values and ligand–protein interactions details are provided in Table 3 and Table S9, respectively. Analysis of the docking results elucidates that compound D1 exhibited a high binding affinity (−9.9 kcal.mol⁻¹), forming stable interactions with the key residues of the active site of EBOV-GP. The interactions formed by compound D1 were three hydrogen bonds (with GLY B:45, TYR B:16, and VAL A:176), one Pi-Sulfur interaction with MET B:47, and three Alkyl and Pi-Alkyl interaction involving (ALA A:74), (LEU B:14), and (VAL A:39) residues. Compound D2 exhibited stability through three hydrogen bonds (two conventional with (PRO A: 160) and (TYR B:16), one carbon hydrogen bound with (LEU B:46), one pi–sulfur interaction with (MET B:47), one pi-pi Stacked interaction with (TYP B: 16), and seven alkyl and pi-alkyl interactions. Similarly, compound D3 formed interactions with (ARG A:37) and (VAL A:176) via conventional hydrogen bonds and with (GLY B:45) via carbon hydrogen bond. Furthermore, it engages in Pi-Donor hydrogen bonds with (TYR B:16), contributing to its structural stability. Additionally, D3 exhibits hydrophobic interactions with (ALA A:74), (LYS A:163), (LYS A:164), (VAL A:39), and (LEU B:14), further enhancing its binding affinity and molecular interactions. Compound D4 was stabilized by the formation of one carbon hydrogen bond with (LYS A:163) and eight hydrophobic interactions with (ALA A:74), (VAL A:39), (LEU A:41), (LEU B:14), (MET B:47) and (PHE A:167). Compound D5 formed one conventional hydrogen bond with (THR B:18) and seven hydrophobic interactions with (ALA A:74), (LEU B:14), (LYS A:163), (MET B:47), (VAL A:39), (THR B: 18), and (TYR B:16).

Close

Fig. 5

Chemical structures of the top five newly designed compounds.

Export to PPT

Close

Fig. 6

2D and 3D representation of the interactions of ligands D1-D5 within EBOV-GP binding pocket.

Export to PPT

3.6 Molecular dynamics simulations

Molecular dynamics simulations were employed to assess the stability of interactions between proteins and ligands, as well as to quantify the binding affinity of ligands. MD simulations were run for the top five compounds (D1-D5). The study included the computation of diverse parameters (Table 4) such as RMSD, RMSF, Rg, and SASA metrics. The best three stable systems (D2, D3, and D5) that showed good stability with Apo form of the EBOV-GP and the EBOV-GP-lead compound complex are presented in the Table 4.

The RMSD serves as a measure to quantify the degree of deviation in a protein's backbone from its initial to final conformation during a simulation, providing valuable insights into its structural stability. A protein with consistent stability will exhibit minimal backbone deviation, while a less stable protein will show greater variation. Analysis of the MD results, as depicted in both the Fig. 7 and Table 4, reveals variations in the RMSD values for the Apo form of EBOV-GP and its complexes with different compounds. Interpreting these findings, the compound with the lowest average RMSD, namely compound D2, suggests a higher degree of stability compared to the complex protein-lead compound. Similarly, compound D3 and compound D5 also exhibit greater stability, as evidenced by their lower RMSD values. This implies that these compounds potentially contribute to maintaining the structural integrity of the protein during the simulation.

Close

Fig. 7

Time-dependent RMSD of c-α backbone of EBOV-GP Apo, lead compound, compound D2, compound D3, and compound D5.

Export to PPT

To further investigate the ligand stability within the binding site of EBOV-GP, the RMSD values of the ligands were monitored (Fig. 8). Monitoring these trajectories would provide valuable information about the conformation and orientation of the simulated ligands within the binding site of EBOV-GP throughout the simulation period. The RMSD plot shows that compound D3 exhibits the most consistent and stable behavior, showing the least variation throughout the 100 ns MD simulation. In contrast, compounds D2 and D5 show higher RMSD values than the lead compound, indicating greater mobility within the binding site before reaching a semblance of stability later at the end of the simulation.

Close

Fig. 8

RMSD of the ligands (reference compound 78 and designed compounds D2, D3 and D4).

Export to PPT

Analyzing the RMSF is essential for identifying flexible areas within protein–ligand complexes. Higher RMSF values in regions with lower structural organization, such as loops and turns, contrast with lower fluctuations in more ordered segments like alpha helices and beta sheets. It is crucial to note that the RMSF assessment covered C-alpha atoms across the entire simulation for both the complexes and the EBOV-GP Apo form. In evaluating the flexibility of the top five compounds binding to the EBOV-GP, RMSF values from the lead compound complex served as a reference baseline. Remarkably, compound D3 demonstrated the lowest RMSF (Table 4 and Fig. 9), corroborating the findings observed in the RMSD analysis. This underscores the consistent behavior of compound D3, not only in terms of structural stability (as indicated by the RMSD) but also in terms of exhibiting minimal fluctuations in flexible regions, further highlighting its potential as a stable and structurally sound ligand for EBOV-GP.

Close

Fig. 9

Plot of RMSF for c-α atoms of EBOV-GP Apo, lead compound, compound D2, compound D3, and compound D5.

Export to PPT

Radius of gyration, serving as an indicator of protein structure compactness, is employed to evaluate changes induced by ligand binding through a comparison of Rg values before and after binding. A lower radius of gyration signifies a more compact protein structure. Fig. 10 presents the Rg plot for the Apo state, lead compound, compound D2, compound D3, and compound D5 unveiled average Rg values of 2.268, 2.293, 2.312, 2.291, and 2.277 nm, respectively. Particularly noteworthy is that compound D5 exhibited Rg values slightly lower than those of the EBOV-GP bound to the lead compound. This suggests that compound D5 may lead to a protein structure that is more densely packed compared to the lead compound, implying potential implications for the binding process and overall stability of the protein–ligand complex.

Close

Fig. 10

Plot of Rg vs. time for EBOV-GP Apo, lead compound, compound D2, compound D3, and compound D5.

Export to PPT

The SASA parameter refers to the surface area of a protein that engages with solvent molecules, measuring the fraction of a protein's surface available to such molecules.

The SASA values for the EBOV-GP Apo state, lead compound, compound D2, compound D3, and compound D5 are 198.266, 205.813, 206.417, 203.612, and 204.020 nm², respectively. The Apo structure inherently represents solvent accessibility, and the marginal increase in SASA with the lead compound indicates alterations upon ligand binding. Notably, compound D2, as showed in Fig. 11, exhibits higher SASA values, suggesting modifications in the protein's surface accessibility, likely attributed to the binding of distinct ligands. These SASA variations offer insights into the protein's dynamic nature, highlighting structural adaptations and potential functional implications associated with specific ligand interactions.

Close

Fig. 11

The SASA profile of EBOV-GP Apo, lead compound, compound D2, compound D3, and compound D5.

Export to PPT

3.7 ADMET

In drug discovery and development, the ADMET emerges as a pivotal phase and plays a vital role in circumventing potential complications during clinical trials. A successful drug candidate must not only exhibit efficacy against the targeted therapeutic endpoint but also demonstrate favorable drug-like qualities and pharmacokinetic characteristics. In the present study, in silico ADMET analysis results of our newly compounds are summarized in Table 5. All the five compounds exhibit a reassuring absence of risks of tumorigenicity, irritancy, mutagenicity, reproductive effects, carcinogenicity, or hepatotoxicity, underscoring their potential safety and efficacy in the drug development pipeline. The human intestinal absorption (HIA) of an oral drug is a critical determinant of its apparent efficacy. The scale of the five compounds ranges from excellent (0–0.3) to poor (0.7–1.0), as it is represented in red in Table 5. Caco-2 cell permeability, expressed as a log cm/s value, is a critical index for evaluating the potential in vivo drug permeability of a compound, with a predicted value >-5.15 log cm/s indicating adequate permeability across intestinal cell membranes via various transport processes. With the exception of compounds D3 and D5, all compounds described have shown good predicted permeability (Barrett et al., 2022). In addition, solubility is a critical factor that influences the ADME of a drug candidate. The solubility is often expressed as the logarithm of its solubility in water (logS). A more negative logS value indicates lower solubility, while a less negative (or more positive) logS value indicates higher solubility. Compared to the lead compound, all the new compounds have higher solubility than the starting lead compound (−5.671), as indicated by their less negative logS values. This suggests an improvement in solubility, which is beneficial for drug development since higher solubility typically correlates with better bioavailability. D5, with the lowest permeability, had the highest solubility value (−4.409) as the permeability is inversely correlated with the solubility (Barrett et al., 2022). For drugs targeting the central nervous system (CNS), traversing the blood–brain barrier (BBB) is imperative to reach their molecular targets, while drugs with peripheral targets may require minimal or no BBB penetration to avoid CNS side effects. The evaluation criteria for BBB penetration are categorized as excellent (0–0.3), medium (0.3–0.7), and poor (0.7–1.0). As detailed in the table below all the proposed structures had excellent to medium values of BBB. Another essential parameter in ADMET analysis is the half-life of a drug (T_1/2), a composite measure of both clearance and volume of distribution. It is a better assessment with reliable estimates of these two properties. Excellent categorization was attributed to all molecules as they are in the range 0–0.3. Additionally, these drug-like molecules are easy to synthesis according to their excellent Synthetic Accessibility scores (SA score ≤ 6). This score was estimated based on a combination of fragment contributions and a complexity penalty (Ertl and Schuffenhauer, 2009). Furthermore, neither inhibitory nor substrate interactions with the human cytochrome P450 family CYP450- 2C9 were observed, suggesting a favorable safety profile with minimal impact on drug metabolism mediated by CYP450-2C9, which is a positive characteristic in drug development. Regarding toxicity, all compounds pose no risk except compound D5, which may have an irritant effect.