SMILES-based QSAR virtual screening to identify potential therapeutics for COVID-19 by targeting 3CL^pro and RdRp viral proteins

Abstract

The COVID-19 pandemic has prompted the medical systems of many countries to develop effective treatments to combat the high rate of infection and death caused by the disease. Within the array of proteins found in SARS-CoV-2, the 3 chymotrypsin-like protease (3CL^pro) holds significance as it plays a crucial role in cleaving polyprotein peptides into distinct functional nonstructural proteins. Meanwhile, RNA-dependent RNA polymerase (RdRp) takes center stage as the key enzyme tasked with replicating the viral genomic RNA within host cells. These proteins, 3CL^pro and RdRp, are deemed optimal subjects for QSAR modeling due to their pivotal functions in the viral lifecycle. In this study, SMILES-based QSAR classification models were developed for a dataset of 2377 compounds that were defined as either active or inactive against 3CLpro and RdRp. Pharmacophore (PH4) and QSAR modeling were used for the virtual screening on 60.2 million compounds including ZINC, ChEMBL, Molport, and MCULE databases to identify new potent inhibitors against 3CL^pro and RdRp. Then, a filter was established based on typical molecular characteristics to identify drug-like molecules. The molecular docking was also performed to evaluate the binding affinity of 156 AND 51 potential inhibitors to 3CL^pro and RdRp, respectively. Among the 15 hits identified based on molecular docking scores, M3, N2, and N4 were identified as promising inhibitors due to their good synthetic accessibility scores (3.07, 3.11, and 3.29 out of 10 for M3, N2, and N4 respectively). These compounds contain amine functional groups, which are known for their crucial role in the binding interactions between drugs and their targets. Consequently, these hits have been chosen for further biological assay studies to validate their activity. They may represent novel 3CL^pro and RdRp inhibitors possessing drug-like properties suitable for COVID-19 therapy.

Introduction

COVID-19, caused by the SARS-CoV-2 virus, has affected over 768 million people worldwide and resulted in more than 6 million deaths [1]. The pandemic has highlighted the urgent need for effective therapeutics to combat emerging viral diseases, given its devastating impact on global health and the economy. While the rapid development of vaccines mitigated the severity of the outbreak, a continued demand persists for antiviral treatments to manage infections, reduce transmission, and address new viral variants. As a result, drug discovery remains a central focus of ongoing research efforts aimed at combating both current and future viral outbreaks.

One promising approach in drug discovery is targeting specific viral proteins essential for viral replication and conserved across different strains of coronaviruses [2]. Notable targets include the 3-chymotrypsin-like protease (3CL^pro) and RNA-dependent RNA polymerase (RdRp), both of which play crucial roles in the replication cycle of SARS-CoV-2 [3, 4]. The 3-chymotrypsin-like protease (3CLpro) is a crucial enzyme in the SARS-CoV-2 life cycle, responsible for cleaving the viral polyprotein into functional components essential for replication. Its unique catalytic Cys-His dyad and conserved active site, capable of accommodating multiple substrates, make it an ideal target for antiviral drug development [5]. Similarly, the RNA-dependent RNA polymerase (RdRp), another key viral enzyme, facilitates viral RNA synthesis and is highly conserved across coronaviruses, making it a pivotal target for therapeutics [6]. Together, 3CLpro and RdRp play essential roles in viral replication, positioning them as prime targets for developing COVID-19 treatments. Targeting these proteins has shown promise in identifying therapeutic candidates, but conventional drug discovery methods face significant limitations, including high costs, time inefficiency, and low success rates [7, 8]. Consequently, more rapid and cost-effective approaches are needed to identify potential antiviral compounds while ensuring their efficacy and safety [2].

Several anti-RNA polymerase drugs currently available, such as Ribavirin [9], Galidesivir [10], Remdesivir [11], and Tenofovir [12], have been approved for treating various viral infections. These drugs are now being evaluated for their effectiveness against SARS-CoV-2 RNA-dependent RNA polymerase (RdRp). Regarding the 3CLpro target, numerous studies and ongoing clinical trials have highlighted drugs like Lopinavir [13], Darunavir [14], Ritonavir [15], Ganovo [16], and Cobicistat [17]. Among these, the combination of Ritonavir/Lopinavir (LPV) is frequently tested in clinical trials for COVID-19 treatment. While there is some evidence suggesting LPV’s potential efficacy, its significant side effects are a major concern [18, 19]. Moreover, these findings underscore the importance of RdRp and 3CLpro as crucial targets for drug development against SARS-CoV-2, with inhibiting their activity emerging as a promising therapeutic strategy.

To address these challenges, quantitative structure–activity relationship (QSAR) machine learning models have emerged as a promising tool for accelerating drug discovery. QSAR models can predict the activity of compounds against specific targets based on their molecular properties, enabling the rapid screening of large compound libraries. These models have been successfully applied in identifying potential therapeutics for various diseases, including cancer and Alzheimer’s [20,21,22].

In QSAR classification, machine learning algorithms classify compounds based on their chemical structure and predicted activity. The goal is to build a model that can accurately predict a compound’s activity against a target by analyzing its structural features, such as molecular weight and hydrophobicity. The model is trained on a dataset of compounds with known activity, which is divided into training and test sets. The performance of the model is evaluated using metrics like sensitivity, specificity, and accuracy, ultimately leading to a binary classification (active or inactive) for each compound in the test set [23, 24].

A particularly advantageous approach is using CORAL (Consensus Modeling for Assessing Chemical Toxicity) and SMILES-based QSAR models. These models can handle large and diverse chemical datasets efficiently, thanks to their use of the Simplified Molecular Input Line Entry System (SMILES) notation. This molecular representation captures key structural features while reducing computational costs, making these models well-suited for virtual screening and prioritizing compounds for experimental testing. The CORAL model, which integrates multiple QSAR models and descriptors to generate consensus predictions, further enhances accuracy and robustness [25,26,27].

In this study, we leverage a dataset compiled by Ivanov et al. [28], which includes compounds tested for activity against the viral targets 3CL^pro and RdRp. By applying QSAR machine learning models, we aim to identify additional compounds with the potential to serve as therapeutics for COVID-19 and related viral infections. An overview of the steps of the present study is shown in Fig. 1.

Fig. 1

Process of molecular modeling in the present study

Materials and methods

Data collection

2377 molecules were selected from an article published by Julian Ivanov and colleagues in 2020 to investigate their inhibitory potency on the COVID-19 virus proteases [28]. The Simplified Molecular Input Line Entry System (SMILES) strings and IC50 values used for CORAL input were directly retrieved from the supplementary files of this article. Initially, the CORAL software checks for duplicated chemicals, incorrect SMILES, and inconsistencies in activity data notation based on SMILES. The SMILES format of the compounds was presented in Table S1 and S2. These molecules include 1168 molecules for 3CL^pro, consisting of 468 active and 700 inactive molecules, and 1209 molecules for RdRp, consisting of 464 active and 745 inactive molecules. Compounds with an IC₅₀ of ≤ 10µM were classified as active compounds, and compounds with an IC₅₀ of ≥ 10µM were classified as inactive compounds. In this study, “semi-correlation” were constructed for ten distinct divisions, as opposed to traditional correlations (Fig. 2)[29, 30]. The splits’ identities were also calculated and are displayed in Table 1. It’s important to mention that there are no pairs of splits with an identity exceeding 35%.

Fig. 2

Visualization of the general concepts of the traditional correlation and semi-correlation

Table 1 Percentages of identity for random splits

CORAL method

The chemical elements in the molecular structure were encoded as symbols for cycles and branching using SMILES attributes, and CORAL software was utilized to build models based on this representation [31]. The CORAL software, which can be downloaded for free from http://www.insilico.eu/coral, is a computational tool that utilizes Monte Carlo methods to develop regression and classification models based on SMILES descriptors and their corresponding endpoints (i.e. active or inactive). To build the models, the compounds were randomly divided into four datasets: a training set (30%), an invisible training set (30%), a calibration set (20%), and a validation set (20%). To create the QSAR model, data from the training set was utilized [30]. The training set is the foundation for constructing the model or “builder”. Monte Carlo optimization is being employed to adjust correlation weights for molecular features derived from SMILES associated with this particular set. The invisible training set is the “inspector” of the model. Computing descriptors for SMILES within this set should either confirm or reject the model's appropriateness for substances not directly engaged in the optimization procedure. The calibration set should identify the onset of overfitting. Based on computational experiments, it's evident that optimization enhances the correlation between descriptors and an endpoint for both the training and invisible training sets. However, as the optimization progresses through more epochs, there's a gradual decrease in the correlation coefficient between descriptors and the endpoint for the calibration set [32]. The validation set was used to test the predictability of the QSAR model. The CORAL software utilized identical algorithms to compute the classification models. The regression models were established on genuine correlations, while the classification models relied on pseudo correlations. The dataset for the classification models was allocated either a value of 1, indicating “active,” or a value of 0, indicating “inactive.”

Optimal descriptor

Molecular descriptors, which are mathematical entities, encode the chemical and physical properties of molecules. Choosing an appropriate set of descriptors is crucial as it determines the accuracy of predictive models. 2D representations based on molecular graphs contain around 70–80% of the information in QSAR models. Fragment-based QSAR calculates descriptors for molecular fragments generated by predefined rules. 2D QSAR is simpler and less time-consuming than other methods. Useful fragments can be identified and mapped onto molecules to suggest improvements.

CORAL software offers three distinct descriptor of correlation weights (DCW) types: graph-based, SMILES-based, and hybrid [33]. The graph-based and SMILES-based DCWs are calculated solely based on their respective input types, while the hybrid DCW is generated using both graph and SMILES inputs. To compute the optimal descriptors using SMILES, the following formula is utilized:

$${}^{SMILES}DCW= \sum_{k=1}^{N}CW\left({S}_{k}\right)+\sum_{k=1}^{N-1}CW\left({SS}_{k}\right)+ \sum_{k=1}^{N-2}CW\left({SSS}_{k}\right)+ \sum_{k}CW\left(NOS{P}_{k}\right)+\sum_{k}CW\left(HAL{O}_{k}\right)+ \sum_{k}CW\left(BON{D}_{k}\right)+\sum_{k}CW\left(HAR{D}_{k}\right)$$

(1)

Table 2 contains a detailed description of the SMILES attributes invariants used in Eq. (1).

Table 2 The detailed description of SMILES attributes

In the context of SMILES fragments, the threshold is a parameter used to distinguish between rare and non-rare fragments. If the number of SMILES containing a particular fragment in the training set is less than the threshold, it is considered rare; otherwise, it is classified as non-rare. During the Monte Carlo optimization process, one epoch refers to a complete cycle of modifying all correlation weights. The correlation weights, which determine the maximum value of the correlation coefficient between the endpoint (0,1) and DCW (Threshold, Nepoch), are obtained numerically using the Monte Carlo method.

The target function of the optimization is to find the correlation coefficient between the optimal descriptor and an endpoint. This is achieved through the CORAL model, which establishes a linear relationship between a predicted endpoint Y and a descriptor of correlation weights (DCW). The DCW is represented by the following mathematical equation:

$$Y={C}_{0}+{C}_{1}{\times} DCW\left({T}^{*},{N}^{*}\right)$$

(2)

Using the least squares method, we can obtain two regression coefficients, C0 and C1, which are used to create an optimal DCW model based on the dataset. During the Monte Carlo optimization process, the parameters T* and N* are used, where T* represents the threshold value and N* represents the number of epochs. In DCW modeling, T refers to the threshold number of active attributes. For example, if T is set to four, any attributes that appear in less than four molecules are considered inactive, and their correlation weight will be zero [34,35,36].

Statistical criteria

For the binary classification model where “active” is represented by 1 and “inactive” by 0, several statistical measures are used to evaluate the model's performance. These measures include the Matthews correlation coefficient (MCC), sensitivity, specificity, and accuracy, which are calculated as follows:

$$Sensitivity=\frac{TP}{TP+FN}$$

(3)

$$Specificity=\frac{TN}{TN+FP}$$

(4)

$$Accuracy=\frac{TP+TN}{TP+FP+FN+TN}$$

(5)

$$MCC=\frac{TP{\times} TN-FP{\times} FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}$$

(6)

A confusion matrix is a table used to assess the performance of a classification model. It includes the number of true positives (TP), true negatives (TN), false positives (FP), and false negatives (FN) predictions made by the model [30].

The performance measures used in evaluating a classification model are important indicators of its quality. One such measure is the Matthews correlation coefficient (MCC), which is commonly used for binary classification problems [37]. The MCC is similar to the traditional correlation coefficient but is designed specifically for this case. The MCC coefficient, applied in machine learning as a balanced measure of the quality of binary classifications, proving beneficial even when the classes vary significantly in size. In practical terms, a model is considered good if the MCC 1, or generally, the MCC should be greater than 0.6 for it to be effective. Another important measure is sensitivity, which assesses the model’s ability to correctly identify positive observations, such as active compounds. Specificity, on the other hand, evaluates the model’s ability to accurately identify negative observations, such as inactive compounds. Finally, accuracy is a measure of the overall performance of the model, taking into account its ability to predict both positive and negative observations. These measures are crucial in determining the effectiveness of a classification model and ensuring its ability to make accurate predictions.

Applicability domain

The OECD QSAR validation principles require QSAR models to be used within their applicability domain (AD), which refers to the space or knowledge used to develop the model and make predictions for new compounds. Regarding CORAL models, the domain of applicability is determined by analyzing the statistical defects of SMILES, which are calculated based on the distribution of available data in the training, invisible training, calibration, and validation sets [38]. To define the domain of applicability, the defect of the SMILES attribute is measured as the difference between the probability of the attribute in the training set and its probability in the calibration set. The total SMILES-defect is obtained by summing up the defects of all attributes. If the SMILES-defect of a particular SMILES is less than double the average defect of compounds in the training set, it is considered to be within the domain of applicability; otherwise, it falls outside the domain of applicability.

$${Defect}_{{A}_{k}}= \frac{\left|{P}_{TRN}(\right.{A}_{k})-\left.{P}_{CAL}\left({A}_{k}\right)\right|}{{N}_{TRN}\left({A}_{k}\right)+ {N}_{CAL}\left({A}_{k}\right)}, if \,{N}_{TRN}\left({A}_{k}\right)>0 \,{Defect}_{{A}_{k}}=1, if \,{N}_{TRN}\left({A}_{k}\right)=0$$

(7)

The probabilities of attribute A in the training and calibration sets are denoted by ${P}_{TRN}\left({A}_{k}\right)$ and ${P}_{CAL}\left({A}_{k}\right)$, respectively. On the other hand, the frequencies of A in the training and calibration sets are represented by ${N}_{TRN}\left({A}_{k}\right)$ and ${N}_{CAL}\left({A}_{k}\right)$, respectively.

In the SMILES notation, the statistical defect (D) is the total sum of all statistical defects of all attributes.

$${Defect}_{Molecule}=\sum_{k=1}^{NA}{Defect}_{{A}_{k} }$$

(8)

In CORAL, the number of active SMILES attributes in a compound is denoted by NA. If a compound is outside the domain of applicability, it is considered an outlier. The detection of outliers in CORAL is based on the condition expressed in inequality 9.

$${Defect}_{Molecule} >2{\times} \overline{{Defect }_{TRN}}$$

(9)

The mean value of statistical defects calculated for the training dataset is referred to as ${Defect}_{TRN} D.$

Virtual screening workflow

Virtual screening is a computational technique used in drug discovery and development to identify potential drug candidates through in-silico (computer-based) screening of large libraries of small molecules. The process involves the use of computer algorithms and molecular modeling techniques to predict the potential binding affinity of small molecules with target proteins, and thus, identify molecules with a high likelihood of being effective drug candidates [39]. The Pharmit webserver was used to conduct in silico virtual screening on four databases (ChEMBL, ZINC, MCULE, and MolPort) in building pharmacophore (PH4) model and the identification of drug-like molecules stages [40]. Pharmit is an online interactive server that allows for the virtual screening of various compound databases using pharmacophore models and molecular shapes. The reference article identified several molecules with in vitro inhibitory effects on 3CL^pro, among which 1H-Indole-2-carboxylic acid, 5-fluoro-, 1H-benzotriazol-1-yl ester (A1) with an IC₅₀ of 0.013 was considered the most active compound. Similarly, 3-[Isopropyl(trans-4-methylcyclohexylcarbonyl)amino]-5 phenylthiophene -2 carboxylic acid (A2) was considered the most active compound with in vitro inhibitory effect on RdRp, with an IC₅₀ of 0.009 [28]. To perform the structure-based pharmacophore screening, we utilized these two most active compounds A1 and A2 in the Pharmit server. The QSAR models were selected to conduct virtual screening. Pharmit filters compounds during hit screening based on drug-like properties including number of rotatable bonds, molecular weight, logP, topological polar surface area, number of HBAs, number of HBDs, and number of aromatic groups. The specified ranges for these properties are: MW ≤ 500, rotatable bonds ≤ 10, logP ≤ 5, PSA ≤ 140 Å2, aromatic groups ≤ 5, 2 ≤ HBA ≤ 7, and 2 ≤ HBD ≤ 7. The binding modes of inhibitors, along with the critical molecular interactions inside 3CL^pro and RdRp active sites with PDB ID 6Y2F and 6NUR, were investigated using the Smina [41] molecular docking package. A purification process was performed by removing all heteroatoms and solvent molecules from the structure. Polar hydrogens were then added to the PDB files. All ligands for docking were sketched using Discovery Studio 2020, and assigned gasteiger charges and energy optimization of ligands using the steepest descent algorithm carried out by Open Babel [42]. The details of the molecular docking algorithm in Smina have been explained in previous studies [36, 43]. Visualization and interaction analyses were performed using the Discovery Studio 2020 viewer.

Results and discussion

QSAR models

For the 3CL^pro enzyme, a collection of 1168 molecules were used, with 468 molecules known to be active against the enzyme and 700 molecules known to be inactive. These molecules were split into four different groups: a training set of 339 molecules, an invisible training set of 368 molecules, a calibration set of 230 molecules, and a validation set of 231 molecules. Similarly, for the RdRp enzyme, a group of 1209 molecules were used, with 464 known to be active and 745 known to be inactive. These molecules were also divided into four groups: a training set of 358 molecules, an invisible training set of 386 molecules, a calibration set of 225 molecules, and a validation set of 240 molecules. To begin with, the process of parameter optimization was carried out in order to determine the most suitable values for the threshold and number of epochs (N_epoch). This was done to ensure that the model would produce accurate predictions while minimizing the risk of overfitting.

To analyze the data set for proteases 3CL^pro and RdRp, it was divided into ten separate parts referred to as “splits” (split 1, split 2, and so on up to split 10). The analysis was conducted using a N_epoch value of 30 and threshold values ranging from 1 to 3. For all ten splits of both proteases, the Monte Carlo optimization used a preferred threshold value (T*) of 1 and a preferred number of epochs (N*) of 3. The models for each of the ten splits (splits 1 through 10) were obtained using this methodology. Tables 3 and 4 displays the statistical properties of the binary classifications for 3CL^pro and RdRp, with the activity being determined using the following formula:

Table 3 the statistical properties of the inhibitor activity classification models of 3CL^Pro that were generated through the Monte Carlo method optimization for ten random splits

Table 4 The statistical properties of the inhibitor activity classification models of RdRp that were generated through the Monte Carlo method optimization for ten random splits

$$Class=\left\{\begin{array}{c}Inhibitory \,activity \le 10\mu M, class=1(active)\\ Inhibitory \,activity >10\mu M, class=0 (inactive)\end{array}\right.$$

(10)

The majority of the models in the study achieved an accuracy, sensitivity, and specificity greater than 90% for the training, invisible training, calibration, and validation sets of 3CL^pro and RdRp across splits 1–10, indicating their ability to predict the activity of the viral enzymes.

It’s important to highlight that the classification model relies on unique semi-correlations, making it unsuitable for using of some traditional criteria. For traditional correlation, an r² value around 0.4 suggests a weak regression model. But when it comes to semi-correlation with a similar r² value (0.4), having specificity, sensitivity, accuracy, and MCC at that level could be seen as good or even outstanding. This means most items in the chart are accurately classified, usually with just one false positive and one false negative. The CORAL software operates under the fundamental assumption that a well-performing model on a calibration set should also perform well on an external validation set. In this case, the statistical parameters of the models for the two proteases are deemed good and acceptable. Furthermore, for split #7 (bold in the Table 3) of 3CL^Pro and split #4 (bold in Table 4) of RdRp, the MCC values for the validation sets are among the top-performing models, with 0.9478 and 0.8890, respectively.

By running Monte Carlo optimization multiple times, one can obtain three groups of SMILES attributes: (i) attributes or promoters that only have positive correlation weights, which promote the activity of compounds; (ii) attributes or promoters that only have negative correlation weights, which promote the inactivity of compounds; and (iii) attributes or promoters that have both positive and negative correlation weights in multiple Monte Carlo runs, and their role is not yet understood. This method allows for a mechanistic interpretation of the model. Tables 5 and 6 provides a list of potential promoters that be linked to both increased and decreased activities for both proteases. The findings from Table 5 revealed that certain structural features, such as aliphatic oxygen and double bond, nitrogen and double bond, aliphatic oxygen and double bond with branching, two successive aliphatic carbon, aliphatic carbon and double bond with branching, double bond, presence of at least three rings, carbon and double bond with branching, aliphatic carbon and aliphatic nitrogen with branching, aliphatic carbon and aliphatic oxygen with branching, molecule containing nitrogen and oxygen with at least one ring and branching, were identified as significant factors for increased inhibitory activity against 3CL^Pro. On the other hand, decreased activity against 3CL^Pro was associated with certain structural features, such as nitrogen and oxygen with double bond, oxygen and double bond, three successive aliphatic carbon, aliphatic carbon and double bond with at least three rings, hydrogen and stereo specific bonds.

Table 5 A collection of SMILES attributes that can be interpreted as promoters of activity or inactivity of compounds against 3CL^Pro

Table 6 A collection of SMILES attributes that can be interpreted as promoters of activity or inactivity of compounds against RdRp

Table 6 presented our findings on the influential features that affect inhibitory activity against RdRp. We observed that negative charge, nitrogen and double bond, aliphatic carbon with branching, double bond, nitrogen and oxygen, double bond, aliphatic carbon and aliphatic nitrogen, aliphatic carbon and branching, presence of at least two rings, the highest number of sulfur equal to zero, oxygen and double bond, and aromatic carbon in the first ring had a positive impact on inhibitory activity. On the other hand, aliphatic carbon with branching and aliphatic oxygen, aliphatic oxygen and two branching, aliphatic carbon and branching, aliphatic nitrogen with branching and aliphatic carbon, branching with aliphatic carbon and stereo specific bonds, successive two aliphatic carbon and aliphatic nitrogen, aliphatic carbon and double bond, aromatic nitrogen in the second ring, the highest number of oxygens equal to one, two successive aliphatic carbons and ring, carbon and double bond with branching, and presence of three successive aliphatic carbon were found to decrease inhibitory activity against RdRp.

Virtual screening analysis

The Pharmit server was utilized to propose pharmacophore characteristics by emphasizing on the key residues involved in active site interactions of 3CL^Pro and RdRp with compound A1 and A2 as the most active compounds, respectively (Fig. 3). Among the databases supported by Pharmit, we opted to screen the ZINC, CHEMBL32, MCULE, and MolPort databases due to the availability of purchasable compounds for virtual screening. The QSAR models were selected for predicting the activity (active or inactive) of compounds. Then, the filter was established using typical molecular characteristics for recognizing drug-like molecules. These features comprise molecular weight, log P (a measurement of lipophilicity), topological polar surface area (a sign of the compound's ability to penetrate cell membranes), the number of rotatable bonds, the number of aromatic groups, the number of hydrogen bond acceptors, and the number of hydrogen bond donors. OpenBabel is utilized to precompute these features. Smina was used in the third screening to confirm that the search for the binding site of 3CL^Pro and RdRp had produced all hits. The compounds were then sorted based on their docking score values, and only those with scores higher than compounds A1 with 3CL^Pro (− 7.22 kcal/mol) and A2 with RdRp(− 7.84 kcal/mol) were included in the list. Figure 3 shows the PH4 models and interaction patterns of two complexes: one between 3CL^Pro and A1, and the other between RdRp and A2.

Fig. 3

The interaction patterns and superpositions of the PH4 models for 3CL^Pro and RdRp with PDB ID (A) 6Y2F with A1 and (B) 6NUR with A2. PH4 features are represented using spheres colored green for hydrophobic, orange for hydrogen bond acceptor, white for hydrogen bond donor, and red for negative ion features

The results of molecular docking indicate that the A1 molecule forms hydrogen bonds with three different amino acids, namely ARG188, GLN192, and THR190. Additionally, this molecule exhibits hydrophobic interactions with amino acids PRO168, MET165, and HIS41. The amino acids GLN189 and VAL186 are also found to be involved in van der Waals interactions with the A1 molecule. The findings of molecular docking analysis between RdRp protease with molecule A2 indicated that GLN724 amino acid forms a hydrogen bond with the carbonyl group's oxygen. Additionally, LEU708 and ARG721 were observed to engage in hydrophobic interactions with the rings and branches, while THR710, GLY712, HIS725, and TYR732 were found to participate in van der Waals interactions. The results are consistent with the pharmacophoric models.

Figure 4 presents the hits obtained from virtual screening procedures across all four databases for each protein. During the final stage of screening, a certain number of compounds were subjected to docking simulations using different proteins, with 156 compounds docked into 3CL^Pro and 51 compounds into RdRp. The structures and molecular docking minimized affinity values of these hits were shown in Table 7.

Fig. 4

Result diagram screening of inhibitors for (A) 3CL^Pro and (B) RdRp from the databases

Table 7 Hits retrieved from the virtual screening alongside their minimized affinity values

Specifically, the hit structures M3 for 3CL^pro, attributed to the presence of rings and oxygen with a double bond, were found to enhance the compound's activity in the QSAR model. Additionally, N2 and N4 for RdRp, characterized by the presence of activity-promoting factors such as nitrogen and a ring in the chain, and the absence of sulfur, are highly compatible with the QSAR models. As a result, they were studied in greater detail. Upon analyzing the interaction between M3 and 3CL^pro, it was observed that there were significant hydrogen bonds formed between the oxygen and nitrogen of M3 with THR25, LEU141, GLY143, SER144, and CYS145. Additionally, two hydrophobic interactions occurred between the two rings of M3 and MET49 and CYS145 (Fig. 5A). According to Fig. 5B, the N2 molecule bound to RdRp forms three hydrogen bonds from N and O atoms with ASN713, GLN724, and TYR732 amino acid residues. Furthermore, the three rings of the molecule engage in hydrophobic interactions with specific amino acid residues such as VAL128, HIS133, LEU207, LEU240, GLY712, and LEU708. The docking results of N5 with RdRp show that there is a hydrogen bond between N5 and TYR732. Hydrophobic interactions were observed between N5 and VAL128, TYR129, HIS133, LEU240, LEU708, and TYR728. Furthermore, ARG132 and ASP465 were involved in attractive charge and salt bridges with the charged part of the molecule (Fig. 5C). The M3, N2, and N4 compounds contain at least an amine functional group, which is known to play a crucial role in drug-target binding interactions. Amine groups can pass through cell membranes when in their non-ionized form, and when ionized, they exhibit favorable solubility in water. Furthermore, the fact that the three compounds are highly soluble in water makes them ideal for various drug production activities. It is worth noting that solubility is a critical feature that affects absorption, especially for discovery projects that aim for oral administration. Additionally, a drug intended for injection must have excellent water solubility to ensure that a small medicinal dose contains enough of the active substance. The three compounds possess a significant number of nitrogens with non-bonding electrons, which suggests that they exhibit hydrophilic properties. As a result, they are logical choices for inhibiting 3CL^pro and RdRp and can be considered as potential candidates.

Fig. 5

The interaction patterns models for (A) 6Y2F with M3, (B) 6NUR with N2 and (C) 6NUR with N5

ADMET properties

ADMET stands for Absorption, Distribution, Metabolism, Excretion, and Toxicity. These are critical factors that affect the pharmacokinetics (how drugs move through the body) and pharmacodynamics (how drugs interact with the body) of a drug. The SwissADME server [44] and the DataWarrior [45] software were used to compute various measures, including the bioavailability score, gastrointestinal absorption, logKp for skin permeation, water solubility and toxicity. The Abbott bioavailability score is used to assess drug-likeness, with a score of 0.55 indicating that the best-predicted hits from virtual screening passed the rule-of-five. Skin permeability, logKp was also calculated and fell within the standard range of − 1 to − 8 for 95% of drugs [46]. The range of water solubility typically reported in drug discovery and development is − 6.5 < logS < 0.5. There are various factors that regulate bioavailability, but ultimately, gastrointestinal absorption is the key determinant [47]. Among the 15 hits identified based on molecular docking scores, M3, N2, and N4 were identified as promising inhibitors due to their good synthetic accessibility scores (3.07, 3.11, and 3.29 out of 10 for M3, N2, and N4 respectively). The reported hits showed high gastrointestinal absorption, and toxicity risk was assessed for mutagenic, tumorigenic, irritant, and reproductive effects. Reproductive toxicity may cause alterations to the male and female reproductive systems, while irritant toxicity can cause reversible damage to the skin or other organs. A majority of the virtual screening hits exhibited satisfactory molecular properties, as shown in Table 8.

Table 8 The projected ADMET characteristics for the identified hits

Conclusion

The COVID-19 pandemic, caused by the novel coronavirus SARS-CoV-2, poses a significant threat to global health. To support the development of effective treatments, we combined ligand-based and structure-based drug design approaches to identify potent inhibitors against SARS-CoV-2. SMILES-based classification models were used to create predictive two-dimensional QSAR models. Our comprehensive virtual screening workflow included PH4 analysis, QSAR modeling, evaluation of drug-like properties, molecular docking, and ADMET testing. The ease of using CORAL software to generate QSAR models proved advantageous for rapidly screening potential compounds, as it reduces the chemical space to a manageable size for further analysis and development. This approach identified 15 potential inhibitors, with M3, N2, and N4 emerging as the most promising candidates due to their favorable synthetic accessibility scores (3.07, 3.11, and 3.29, respectively) and the presence of amine functional groups, which are crucial for effective binding interactions with drug targets. These compounds have been selected for further biological assays to validate their efficacy. This study provides valuable insights into the development of novel inhibitors with potential therapeutic applications for COVID-19. To confirm our computational findings, experimental evaluation of the identified compounds against RdRp and 3CLpro activity using standard enzyme-based or cell-based assays is recommended.