This study proceeded in two phases as illustrated in the graphical depiction in Figure 1. First, a high throughput screen (HTS) measured in Cell-Based and Microorganism Combination System bioassay data of 343,305 compounds retrieved from PubChem was used for the study. The dataset with PubChem AID: 651640 was aimed to identify inhibitors of Dengue Virus by treating BHK-21 with various compounds before being infected with Dengue Virus serotype 2. The DENV2 strain New Guinea C, obtained from the American Type Culture Collection (ATCC) with catalog number VR-1584 was adapted (Che et al., 2009). Several viruses have been successfully studied using the cytopathic effect (CPE) assay to find new antiviral substances (Lin et al., 2023; McCormick et al., 2012).The expected outcome was that compounds increased ATP levels, indicating increased cell viability compared to the positive control (uninfected cells), and the neutral control was considered positive and otherwise negative. The dataset retrieved for this study comprised 5,946 active and 321,638 inactive compounds, a ratio of about 1:50 compounds. To compute molecular descriptors, the actives and inactives of the dataset were downloaded in the Structure Data Format (SDF).

Recent studies have demonstrated the utility of molecular descriptors in predicting compound activity, toxicity, and other properties from chemical structures (Comesana et al., 2022; Trinh et al., 2023). The active and inactive datasets were converted from SDF to Simplified Molecular-Input Line-Entry System (SMILES) formats for easier manipulation for machine learning. PaDEL-descriptor calculator (Yap, 2011) was used to compute 1440 molecular descriptors from the canonical SMILES of the compounds. Before descriptor calculations, compound standardization was performed, involving salt removal and nitro group normalization, to ensure uniformity and accuracy of the descriptor data (Viganò et al., 2024). The dataset used for training and testing the models comprised 4470 active compounds and 16780 inactive compounds, after reducing the ratio between the active and inactive compounds to 1:4 to enhance the computational efficiency.

The significant imbalance between active and inactive compounds was acknowledged, but oversampling techniques to balance the classes were not applied. A variance filter was employed to address the issue of dimensionality and enhance the relevance of features for predictive modeling (Velliangiri et al., 2019). Mean imputation was applied to handle any missing values within the dataset, ensuring a complete set of descriptors for each compound. Standardization of the dataset was done by removing the mean and scaling to unit variance using the Standard Scaler from the scikit-learn library. This transformation ensures that the features are centered around zero and have a standard deviation of one, which is crucial for many machine learning algorithms to perform optimally, as it ensures that the model is not biased towards features with larger numerical ranges. The was computed using the formula 1: Standard Scaler X = X   i − m e a n s t a n d a r d   d e v i a t i o n   (1) Where i represent each value in the feature X.

Five models were built from the training dataset using five machine learning algorithms and the model with superior performance in terms of classification metrics was chosen (Kee et al., 2023; Tougui et al., 2021). A 70%–30% split of the pre-processed data was used for training and testing the models, respectively. The classification models developed included Support Vector Machine (SVM), k-nearest neighbors (k-NN), Gaussian Naïve Bayes (GaussianNB) (Adams et al., 2022), Random Forest classifier (RF), and Logistic Regression (LR) (Khorshid et al., 2021). For the k-nearest neighbors (k-NN) model, k = 3 was used. The Gaussian Naïve Bayes model was implemented using default settings from the scikit-learn library (Pedregosa et al., 2011). The SVM model was optimized with the probability parameter set to True (Sandhu et al., 2022). The Random Forest model was built with a maximum depth of 8 and 100 estimators (J. Adams et al., 2022). Finally, the Logistic Regression model was constructed with a maximum iteration parameter of 1000.

Prior to comparison with other optimized classifiers, each classifier underwent optimization to determine the optimal hyperparameters that yielded the maximum accuracy. The optimized models were evaluated using 10-fold cross-validation. This method splits the training data into kkk groups, trains the model on k-1k-1k-1 folds, and tests the model on the remaining fold to yield a trustworthy estimate of the model’s performance on unseen data (Jung and Hu, 2015; Vabalas et al., 2019). It computes performance metrics such as accuracy, accuracy, precision, recall, and F1 score based on the confusion matrix, which includes true positives (TP), true negatives (TN), false positives (FP), and false negatives (FN) (Tharwat, 2021). These metrics (as shown in Equations 2–6) were used to compare and select the best-performing model for predicting Dengue Virus inhibitors (Adams et al., 2022; Orozco-arias et al., 2020) Accuracy = T P + T N T P + F P + T N   F N (2) Precision = T P T P + F P (3) Recall = T P T P + F N (4) F 1 = 2   x   P r e c i s i o n   ∗   R e c a l l   P r e c i s i o n + R e c a l l (5) Specificity = T N T N + F P (6)

The best-performing model was used to make predictions against 812 compounds from Afrodb, and 1871 EANPDB compounds. Prior to compound predictions, the model’s predictive power was validated on known Dengue Virus inhibitors post-cross-validation. The molecular descriptors of these inhibitors were calculated and preprocessed similarly to the training and test data. This was done to reinforce the credibility of the cross-validation results for the best-performing model based on the metrics. The number of correct predictions made by the model on the submitted inhibitors defines the accuracy of the classification or prediction. A total of 18 known inhibitors curated from literature were submitted to the LR model as a result of their activity against the Dengue Virus shown in different studies.

The Crystal structure of the Dengue 2 Virus nonstructural protein NS2B/NS3 was obtained from the RCSB Protein Data Bank (http://www.rcsb.org/pdb) with PDB ID: 2FOM. Before selecting 2FOM, different suitable structures such as 4M9M, 4M9T, and 4M9I were retrieved, however, 2FOM was chosen based on its low resolution, low R values, and the number of missing residues. The selected protein structure was superimposed with the other suitable structure using PyMOL to measure their root mean square deviation (RMSD). The structure was cleaned using PyMOL (Yuan et al., 2017) to devoid the protein of ions, water molecules, and other structures like ligands before minimization was carried out by employing Groningen Machine for Chemical Simulations (GROMACS) (Abraham et al., 2015). The steepest descent minimization algorithm with a maximum number of 50,000 steps and a minimization step size of 0.01 was used to minimize the protein structure. The three-dimensional structures (.sdf) of the natural compounds predicted by the machine learning model were obtained for the molecular docking stage. Natural products were chosen for this investigation because of their structural and chemical variety, and the therapeutic effects of phytochemicals found in plants.

The molecular docking procedure for the predicted compounds was carried out using AutoDock Vina (Trott and Olson, 2010). All compounds were energy minimized in 200 steps using the Universal force field (UFF) before being translated to the Protein Data Bank partial charge and atom type (.pdbqt) format using Open Babel software (O’Boyle et al., 2011). Visualization of the resultant energy minimized protein structure and the removal of surrounding water molecules before the virtual screening was done in PyMOL v1.5.0.4. The prepared structure was then saved using PyMOL before applying the “make macromolecule” option in PyRx to prepare for the docking of selected hits. The library was screened against the NS2B/NS3 protease using a grid box dimension of center_x = −5.179Å center_y = −9.575Å center_z = 13.756Å size_x = 18.302Å size_y = 19.821Å size_z = 23.788Å.

All hit compounds that contributed binding affinities of at least −8.0 kcal/mol were considered. The output of AutoDock Vina is ranked in a decreasing order of binding affinity using a negative function; a more negative binding affinity is preferred. The best mode for each compound was applied using the root mean squared deviation. The result was then examined with PyMOL to find the optimal-docked ligands. LigPlot + (v1.4.5) was used to analyze the interactions between key residues in the active site of the protein and the docked compounds (Laskowski and Swindells, 2011). The protein-ligand complexes generated in PyMOL were used as input for LigPlot. The resulting output provides a 2D depiction of intermolecular interactions, including hydrophobic interactions and hydrogen bonds.

A pharmacokinetics profile, comprising an assessment of absorption, distribution, metabolism, and excretion (ADME), was applied to a subset of compounds using SwissADME (Daina et al., 2017). Along with Veber’s rule, the ADMET testing also measured five properties: total average molecular weight in g/mol, the number of hydrogen bond donors, hydrogen bond acceptors, rotatable bonds, and partition coefficient—collectively referred to as Lipinski’s rule of five (Devadasu et al., 2018; Jia et al., 2020). Using OSIRIS DataWarrior v06.02.05 (Sander et al., 2015), the toxicity characteristics, including mutagenicity and tumorigenicity, were predicted.

A 100 ns MD simulation was run on each protein-ligand complex and the unbound protein using GROMACS-2020.5 on a Dell EMC high-performance computing cluster at the WACCBIP, University of Ghana, Accra. The CHARMM36 all-atom force field produced the protein and ligand topology (July 2022). Utilizing a cubic box for all simulations, the systems were each solvated with water molecules, neutralized with ions, and energy-minimized to optimize the system. To equilibrate each system, NVT, and NPT ensemble were applied for 100 ps a piece before the production run. The parameters for the production run included 50,000,000 steps which translates into 100 ns and a step size of 0.002 (2 fs). Xmgrace (Turner, 2005) was used to visualize and analyze the root mean square deviation (RMSD), root mean square fluctuation (RMSF), and radius of gyration (Rg) obtained from the MD simulations.

The binding free energies of the protein-ligand complexes and the individual energy contributions of the residues were calculated using the Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) method (Kumari et al., 2014). This is a corroboration technique used to verify the limitations of the existing scoring function utilizing the MD simulation output files (Wang C. et al., 2018). R programming software was utilized to plot the graphs from the MMPBSA computations.

The bioactive dataset obtained from PubChem was an imbalanced set, with approximately one-third of its constituents being active compounds. As displayed in Figure 2, inactive compounds dominated the dataset. Using PaDEL, 1,444 molecular descriptors were generated, which provide a mathematical representation of the compounds for QSAR modeling by converting chemical information about the compounds into numerical values. The dataset underwent a three-dataset splitting strategy to be divided into training, validation, and test sets. The dataset of 21,250 compounds was split into 14,875 training data, 3,187 test data, and 3,188 externally held data. The first set is for training the algorithm, the validation set tunes hyperparameters, and the test set is used to test model performance and predicting ability. The application of a variance filter reduced the number of descriptors from 1,444 to 684, using a variance threshold of 0.1 to filter out descriptors with minimal variance. This step ensured that only the most informative features were retained for subsequent modeling, as low-variance features are often less useful in distinguishing between classes (Velliangiri et al., 2019).

Five machine learning models were developed to predict Dengue Virus inhibitors: k-Nearest Neighbors (k-NN), Gaussian Naïve Bayes (NB), Support Vector Machine (SVM), Random Forest (RF), and Logistic Regression (LR). Each model was evaluated using statistical metrics such as accuracy, precision, recall, and F1 score, as displayed in Table 1 with the confusion matrix shown in Supplementary Table S1.

Among the models, Logistic Regression demonstrated the highest performance across most metrics, including an accuracy of 94%, a precision of 91%, and an F1 score of 0.83. This superior performance made Logistic Regression the most suitable model for predicting Dengue Virus inhibitors in this study. The SVM model followed closely, with an accuracy of 93%, a precision of 94%, and an F1 score of 0.81. Although the k-NN and Random Forest models performed well, they lagged behind the top two models in recall and F1 scores, indicating that they were less effective in identifying all active compounds.

Gaussian Naïve Bayes was the poorest-performing model, with an accuracy of 81% and an F1 score of 0.51. This model’s low precision (0.55) and recall (0.47) indicate that it struggled to balance identifying true positives and minimizing false positives, particularly in the imbalanced dataset where inactive compounds were predominant. Overall, Logistic Regression emerged as the most reliable model for predicting potential inhibitors due to its robust performance across the different metrics.

The study further validated the performance of the developed models by testing 18 known Dengue Virus inhibitors curated from literature. The prediction made by the Logistic Regression model for these compounds is listed in Table 2. The Logistic Regression model correctly predicted 11 of the 18 inhibitors as active, outperforming other models, including SVM and Random Forest.

Several inhibitors, such as Pentoxifylline, Prochlorperazine, and Balapiravir, were correctly classified as active by Logistic Regression, in line with their established mechanisms of action against the Dengue Virus. Notably, inhibitors like Celgosivir and Bortezomib, which inhibit viral replication, were also predicted accurately. This validation process of known inhibitors provided confidence in the model’s predictive capability, suggesting it could effectively generalize to novel compounds with similar mechanisms of action.

Following this validation, the Logistic Regression model was applied to predict 2,683 new compounds, including 812 from the ZINC database and 1,871 from the EANPDB database. Out of these, 933 compounds were predicted to be active, representing a promising pool of potential Dengue Virus inhibitors for further experimental validation.

These predictions highlight the utility of the developed QSAR models, particularly the Logistic Regression model, in identifying novel drug candidates for Dengue Virus inhibition. The robust performance on known inhibitors and newly predicted compounds highlights its potential as a valuable tool in future drug discovery efforts targeting the Dengue Virus. The study’s ability to handle imbalanced data effectively and generate accurate predictions underscores the importance of appropriate descriptor selection and data preprocessing in QSAR modeling.

In this study, the NS2B/NS3 protease was preferentially selected amongst the seven nonstructural proteins of the Dengue Virus as the target structure to corroborate the prediction by the logistic regression model. The NS2B/NS3 protease is an essential enzyme for viral replication and assembly, making it a principal antiviral target for developing therapeutics against the virus (Erbel et al., 2006; Norshidah et al., 2023). There are two potential locations for inhibiting DENV protease: the active site and the blocking attachment of protease (NS3) to its protein cofactor (NS2B). The active site on the NS3 which is the prime target is made up of a conserved catalytic triad like His51-Asp75-Ser135 (Noble et al., 2012; Zamri et al., 2019). A search via the Protein Data Bank repository for a solved structure of the NS2B/NS3 for the Dengue Virus serotype showed IDs such as 4M9T, 2FOM, 4M9M, and 4M9I with resolutions 1.74, 1.50, 1.53, and 2.40 Å respectively and R-value work of 0.215, 0.176, 0.203, and 0.215 respectively. The 2FOM, solved using x-ray diffraction, was selected for this study since it had the lowest resolution and R-value, both of which are a measure of the quality of the structure (Wlodawer et al., 2008). Additionally, the RMSD values after superimposing the selected protein structure to 4M9T, 4M9M, and 4M9I were 0.358, 0.151, and 0.276 Å, respectively, underlying their close structural similarity with the 2FOM. The three-dimensional structure of the 2FOM with a ligand docked in the active site is shown Figure 3. The active site on the NS3 used in this study, besides consisting of the catalytic triad which is pivotal in inhibiting its activity, also consists of residues such as Leu128, Pro132, Ser131, and Tyr161 (Norshidah et al., 2023; Purohit et al., 2022).

A total of 853 compounds from the Logistic Regression prediction were docked into the active site of the target using Autodock. Anhydrophlegmacin showed the highest binding affinity of −9.2 kcal/mol towards the protease among all docked ligands. The docked compounds demonstrated binding affinities between −9.2 and −3.6 kcal/mol. This reinforces the prediction ability of the Logistic Regression model. Applying a threshold of −8.0 kcal, 59 ligands with affinities of −8.0 kcal/mol or better were considered for downstream analysis. A higher cutoff of −7.0 kcal/mol, taken as the standard threshold for a compound to be considered active against a particular target, was employed (Kwofie et al., 2019a). The higher the binding affinity, the stronger the bond between the ligands and the target protein. The protein in the complex with the ligands was visually inspected using PyMOL (Figure 4) to select the best-docked possess. In addition, inhibitors such as Leflunomide and Prednisolone were also incorporated into the docking to act as a control. They demonstrated affinities of −7.1 and −7.0 kcal/mol respectively. Table 3 shows the binding affinities of the top 20 of the selected compounds and the inhibitors.

In continuation of the structure-based molecular docking to further confirm the binding affinities of the predicted compounds, the interactions of the compounds within the predicted binding pocket were determined. The biomolecular interactions between the NS2B/NS3 protease and the compounds were generated using LigPlot. Studies into these interactions are crucial in determining promising lead compounds. Identification of crucial residues in the corresponding targets’ active sites was made possible by characterizing the binding interactions. To find the suitable compound that inhibits the activities of the NS2B/NS3 protease, hydrogen and hydrophobic interactions between the shortlisted compounds and the residues in the active site were elucidated.

For the interactions of the protease, the ligands docked to the active site were observed to interact with the proposed residues such as His51, Ser135, Leu128, Pro132, Ser131, Tyr161, and Asp75 as shown in Table 3 and Supplementary Table S2. Anhydrophlegmacin and anhydrophlegmacin-9,10-quinones_B2 which had the highest binding affinities interacted with similar residues such as His51, Asp75, Gly151, Leu128, Pro132, and Gly153. They interacted with conserved catalytic triad residues Asp75, Ser135, and His51 through hydrogen bonding with bond lengths of 2.57, 3.06, and 2.86 Å respectively. ZINC38628344, which had an affinity of −8.5 kcal/mol with the NS2B/NS3 protease formed hydrogen bond interaction with His51 (2.89 Å), Ser135 (2.68 Å), Asp75 (2.57 Å), Phe130 (3.06 Å), Tyr150 (3.10 Å) and hydrophobic interactions with residues Pro132, Ser131, Leu128, Tyr161, Gly153, Gly151 (Figure 4). The inhibitor Prednisolone interacted via hydrogen bonding with Gly151 (2.90 Å, 2.71 Å), Asp75 (2.95 Å), His51 (3.21 Å), and Gly153 (2.93 Å, 3.16 Å) as shown in Figure 5. In addition, ZINC14441502 formed hydrogen interactions with Gly151 and Ser135 with bond length 2.86 and 2.99 Å respectively; and hydrophobic bonding with Leu128, Gly153, Asn152, Val72, Asp75, His151 and Phe130 (Supplementary Table S2). 39 out of 56 hits docked firmly and interacted with critical residues in the active site, and these were selected for downstream analysis.

Pharmacokinetics controls how medications are absorbed by the body and eventually eliminated (Wang G. Y. et al., 2018). Analyses were conducted on pharmacokinetic features, including gastrointestinal (GI) absorption. Drugs taken orally can enter the bloodstream through a process known as gastrointestinal absorption (GI) (Suenderhauf et al., 2012). “High” compound absorption occurs in the GI tract. To select druglike compounds, Veber’s criteria were also applied, and the selected hits that do not conform to Lipinski’s rule of five (RO5) were eliminated (Ogbodo et al., 2023). Out of the 39 hits, 20 were in violation of the rule (see Supplementary Table S3). Twelve hits also broke one of the RO5s. The remaining 7 hits: 5,7′-physcion-fallacinol, ZINC000095485956, ZINC000085594516, amentoflavone, ZINC000095486111, voucapane-18,19-di-(4-methyl)-benzenesulphonate, ZINC000095485927 showed the least drug-likeness of two RO5 violations (Table 4). Veberr’s rule, with TPSA ≤140 and rotatable bonds ≤10, was used as the main determinant (Veber et al., 2002). 26 out of the selected hits demonstrated 0 violations with the remaining showing only one violation of the rule. The solubility and pharmacological profiles such as GI absorption were also elucidated. Only ZINC000095485927 was predicted to be insoluble (Supplementary Table S3). Four of the 14 hits had a moderate solubility prediction and four had a soluble prediction. However, 19 of the selected hits were predicted to be poorly soluble (Supplementary Table S3). Compounds are considered to have met the GI absorption criteria if it is denoted as ‘High’ suggesting a high propensity of absorption into the intestinal tract for orally administered drugs. 21 and 18 of the selected hits were estimated to be high and low respectively. The mutagenicity and tumorigenicity levels of the hits were also predicted using DataWarrior (Table 4). From the results obtained, 26 out of the 39 hits tested were neither mutagenic nor tumorigenic.

Molecular dynamics simulations were carried out using GROMACS 2020.5 to further elucidate the stability of the predicted lead compounds within the active site of the NS2B/NS3 protein (Mazumder et al., 2017). Understanding the binding mechanisms of the various molecules in the active site is crucial for the design of efficacious drugs. To analyze the dynamic behavior of the unbound proteins and complexes, the root mean square deviation, the radius of gyration, and the root mean square fluctuation were plotted with the use of Xmgrace (Agyapong et al., 2021; Kwofie et al., 2019b; Musyoka et al., 2016). All simulations were carried out for 100 ns.