2.1 Data collection and curation

In the current work, a set of 30 sulfonamide-based bisubstrate analogs were collected from the Binding database (https://www.bindingdb.org/) as SARS-CoV-2 Nsp14 inhibitors, and these molecules target the N7-MTase of SARS-CoV-2 Nsp14 with significant inhibitory activity (Table S1). According to OECD principle 1, “a QSAR should be associated with a defined endpoint” (OECD, 2007), the experimental activity of each compound was expressed as IC₅₀ (μM). The IC₅₀ values are determined based on the same bioassay protocol (based on the radioactive N7-MTase assay) according to previous study (Ahmed-Belkacem et al., 2022). Firstly, the compounds were dissolved in 5 % DMSO, and the SARS-CoV-2 Nsp14 proteins were dissolved in 3 % DMSO. At 30 °C, the compounds with different concentrations were incubated with 50 nM SARS-CoV-2 Nsp14 in the reaction mixture [40 mM Tris-HCl (pH 8.0), 1 mM DTT, 1 mM MgCl₂, 2 μM SAM and 0.1 μM 3H-SAM (Perkin Elmer)] for 30 min. The final enzymatic activity (IC₅₀) was determined by filter-binding assay (Ahmed-Belkacem et al., 2022). For model development, we converted the IC₅₀ (μM) values into molar units (M) and then further converted them into pIC₅₀ (pIC₅₀ = -logIC₅₀) values, which is a customary practice in QSAR modelling (Gramatica et al., 2016; Gramatica, 2020; Li et al., 2022, 2023; Chen et al., 2023). The pIC₅₀ value was positively correlated with the inhibitory activity compared to the IC₅₀ value. We carefully examined all structures in the dataset before calculating the descriptors in order to develop statistically significant QSAR models.

2.2 Descriptor calculation and dataset division

After downloading the structures of the modelling compounds from the Binding database, the structures were optimized by the 3D module of ChemBioDraw Ultra 14.0 software in preparation for the next step of descriptor calculation. Dragon (https://chm.kode-solutions.net/) (Dragon, 2021) and PaDEL descriptors (Yap, 2011) were combined to characterize the structural features of the compounds, including functional group counts, topological indices, ring descriptors, connectivity indices, atomic center fragments, 2D atomic pairs, atomic type E states, molecular properties, etc., to more fully extract the hidden important structural information. To facilitate mechanistic interpretation, only 2D descriptors with clear physicochemical significance were calculated for model building, and constant or near-constant and highly correlated descriptors (correlation coefficient > 0.95) were removed to reduce redundant variable information (Gramatica et al., 2013; Gramatica et al., 2014).

To obtain QSAR models with higher prediction quality, a total of six division techniques in QSARINS 2.2.4 software and Dataset Division GUI software (http://dtclab.webs.com/softwaretools) were used to divide the dataset into training and test sets in a ratio of approximately 3:1. Three division techniques in QSARINS, including (1) ordered by response (ORes), (2) ordered by structure similarity (OStr), and (3) randomly (Rnd) (Gramatica et al., 2013). Dataset Division GUI software includes three division techniques, namely, (1) Kennard Stone based method (KS), (2) Euclidean Distance based method (ED), (3) Activity/Property based method (A/P) (Ambure et al., 2019). The training set was used to build the QSAR model and the test set was used to validate the developed model in terms of its predictive power and to avoid possible biases.

2.3 Descriptor selection and model building

According to OECD principle 2, “a QSAR should be associated with an unambiguous algorithm” (OECD, 2007). Identifying important descriptors from a large number of initial descriptors using different variable selection methods is an important step in QSAR modelling. Due to the small dataset in this study, two different software were used to construct the QSAR models in order to filter a more reliable model. First, the Genetic Algorithm-Variable Subset Selection (GA-VSS) module included in the QSARINS software was used for variable screening, and MLR was used to build the QSAR model. In the descriptor selection process, the leave-one-out cross-validation correlation coefficient ( $Q_{\mathit{LOO}}^2$ ) was used as the fitness function to evaluate the model robustness. The QUIK (Q under influence of K) rule was applied to avoid multicollinearity with a threshold of 0.05 (Todeschini et al., 1999). Secondly, the three dataset division results generated by Dataset Division GUI software were further used by the genetic algorithm (GA) in the Small Dataset Modeler software (https://sites.google.com/view/smalldatasetmodelling) as a variable selection method to build the MLR model and its corresponding PLS model (GA-PLS model using the same descriptors selected in the GA-MLR analysis) (Ambure et al., 2019). The final MLR/PLS models were selected based on the MAE_LOO(95%) after running GA-MLR for at least 10 times and the selected model was evaluated using the test set predictions (Roy et al., 2016).

2.4 Model performance evaluation

OECD principle 4 states that “a QSAR should be associated with appropriate measures of goodness-of-fit, robustness and predictivity” (OECD, 2007), and statistical validation is a key step in testing the eligibility of the model. In this study, internationally accepted internal and external validation metrics were used to check the robustness and predictability of the generated models.

In internal validation, we evaluated the robustness and goodness-of-fit of the model mainly through $Q_{\mathit{LOO}}^2$ and the coefficient of determination R² (including adjusted $R_{\mathit{adj}}^2$ ), with larger values indicating better model quality (Tropsha, 2010). As for the external validation, we calculated various statistical parameters to verify the model prediction performance, such as $Q_{F1}^2$ , $Q_{F2}^2$ , $Q_{F3}^2$ , $R_{\mathit{test}}^2$ , $r_{m(test)}^2$ , Δ $r_{m(test)}^2$ , concordance correlation coefficient for the test set (CCC_test), and the root mean square error (RMSE) (Tropsha, 2010; Gramatica and Sangion, 2016). The parameter scales to be satisfied by a better model are as follows (Chirico and Gramatica, 2011, 2012): R² > 0.7, $Q_{\mathit{LOO}}^2$  > 0.6, $Q_{\mathit{Fn}}^2$  > 0.6; CCC_test > 0.85; $r_{m(test)}^2$  > 0.65; Δ $r_{m(test)}^2$  < 0.2; MAE_test and RMSE_test should be as small as possible. For more detailed information on these statistical parameters, please refer to Table S2.

The models generated by QSARINS software were ranked using the multi-criteria decision making (MCDM) module (Gramatica et al., 2013). MCDM_fit was obtained based on the maximum values of R², R_adj², CCC_train and the minimum value of R²-R_adj². MCDM_test was selected based on the maximum values of Q_F1², Q_F2², Q_F3², and CCC_test. Finally, we ranked the generated QSAR models based on MCDM_fit and MCDM_test values. For the models generated by Small Dataset Modeler software (Ambure et al., 2019), the MAE_95% value of the test set of each model was compared, the model with the smallest MAE_95% value was selected.

2.5 Applicability domain (AD)

According to OECD principle 3, “a QSAR should be associated with a defined domain of applicability” (OECD, 2007). For this purpose, the applicability domain (AD) analysis was performed to define the scope of application for a QSAR model. The prediction of a compound can be considered reliable only if it falls within the AD. In this study, we used the leverage method combined with the standardized residual method to define the AD. It can be simply understood that a compound should be identified as a structural outlier if its leverage value h is larger than the threshold value h* (h*=3(p + 1)/n; p: the number of descriptors in the model; n: the number of compounds in the training set); while a compound with predicted residual higher than ± 3 standardized residuals, it should be considered a response outlier (Gramatica et al., 2013). The results of the AD analysis can be visualized by Williams plot, which directly shows the structural and response outliers, but it should be noted that some compounds may not be considered as outliers even if h > h*, because their standardized residual values may be small.

2.6 Screening of untested compounds using developed QSAR models

One of the main purposes of QSAR modelling in this study was to screen compounds having significant inhibitory effects on SARS-CoV-2 Nsp14 from the large database of compounds. We downloaded 263 compounds from the PubChem website (https://pubchem.ncbi.nlm.nih.gov/) using the sulfonamide structural fragment (Fig. 2A) as template molecule 1 to form a large true external set. Notably, the structures of these 263 unknown true external compounds have the same substructures with the 30 compounds used for QSAR modelling. Meanwhile, in order to select compounds with more similar structures to the 30 modelling compounds, a benzene ring was added to template 1 as template 2 (Fig. 2B). Therefore, a total of 42 compounds were extracted from the 263 compounds to form a small true external set. The detailed information of true external set was listed in Table S1. In order to verify the reliability of our chosen QSAR model for the prediction of unknown compounds, certain measures will be chosen to further evaluate the prediction results, which will help us to prioritize potential compounds in expensive in vitro or in vivo experiments. We used “Prediction Reliability Indicator (PRI)” tool (https://dtclab.webs.com/software-tools) to check the prediction confidence of each true external compound (Roy et al., 2018). On the basis of the composite score of each compound, the prediction quality was classified into three groups: “good” (composite score 3), “moderate” (composite score 2) and “bad/unreliable” (composite score 1) (Khan and Roy, 2019; Roy et al., 2018).

Close

Fig. 2

(A) Template 1 for large true external set (263 compounds); (B) Template 2 for small true external set (42 compounds).

Export to PPT

2.7 Molecular docking study

In this paper, we used molecular docking technique to determine the interaction pattern of small molecule ligands with the target protein Nsp14. The crystal structure of Nsp14 protein was retrieved from the Protein Data Bank (PDB) with the code 7R2V (available at: https://www.rcsb.org/structure/7R2V) (Czarna et al., 2022). 7R2V is a protein–ligand complex with S-adenosyl homocysteine (SAH) bound to it. Molecular docking studies were performed using Discovery Studio 4.0 (Accelrys Inc.). Before docking, the ligand SAH bound to the protein was removed, and then the protein preparation was completed by the removal of deletion residues, hydrogenation and docking site generation from the protein. The sites were generated using the “from current selection” program in the receptor-ligand interaction module, referring to the study by Ahmed-Belkacem et al. using the four residues TYR420, PHE506, PHE426, ARG310 to generate active sites (Ahmed-Belkacem et al., 2022). The selected compounds were subjected to ligand preparation through the small molecule module of the Discovery Studio platform to generate a series of ligand conformations. Each of these generated conformations will next be used for molecular docking in the CDOCKER module. The CDOCKER binding energy parameters (kcal/mol) of all the receptor-ligand complexes were examined and the highest scores (most negative) favored the binding of the receptor to the ligand.

2.8 ADMET profiling

After molecular docking study, we tested the screened compounds again to determine their drug-likeness and ADMET properties. We used pkCSM (Pires et al., 2015) (http://biosig.unimelb.edu.au/pkcsm/), SwissADME (http://www.swissadme.ch/) (Daina et al., 2017), and the ADMET and TOPKAT® module embedded in Biovia Discovery Studio software (Kumar et al., 2022). Three different methods were used to predict the ADMET properties of the compounds. The aim of using multiple software is to maximize the assessment of all possible toxicities of the compounds as well as key drug-likeness information, which can increase the chance of reaching anti-SARS-CoV-2 drug candidates.

3.1 QSAR analysis

According to the OECD's QSAR validation guidelines (OECD, 2007), QSAR models that meet various parameter criteria should be selected. In this study, a variety of dataset division methods combined with two computational methods (PLS and MLR) were used to construct nine different classes of models. As shown in Table 1, the nine models were able to explain training set variance (R²) between 0.7150 and 0.7796 (0.6130 ∼ 0.7373 in terms of $Q_{\mathit{LOO}}^2$ ). In addition, $Q_{F1}^2$ ranged between 0.5573 and 0.8070, $Q_{F2}^2$ ranged between 0.5547 and 0.8023, CCC_test ranged between 0.7700 and 0.9073, and MAE_test ranged between 0.2323 and 0.4155. These parameters indicated that all nine models exhibited accepted but with different degrees of predictive ability.

In order to find new potential Nsp14 inhibitors using the best QSAR model and to identify the structural features in the compounds that have important influence on Nsp14 inhibition, we comprehensively compared the validation parameters of the nine models in Table 1. It could be found that the internal and external parameters of Model 1 (Eq.1) were better than those of other models, thus being the best QSAR model. Model 1 (Eq.1) was constructed based on the OStr division method combined with the MLR method. The model contained only two 2D descriptors, namely nCb- and MDEC-12, the correlation between these two descriptors was very weak (Fig. 3A). The typical statistical parameters of Model 1 were R² = 0.7796, $Q_{\mathit{LOO}}^2$  = 0.7373, $R_{\mathit{test}}^2$  = 0.8539, $Q_{\mathit{Fn}}^2$  = 0.8023 ∼ 0.8932 and CCC_test = 0.9073. It is evident that the best model meets the OECD guidelines in terms of robustness, predictivity and adaptability (OECD, 2007). The detailed validation parameters of Model 1 were given in Table 1. The scatter plot (Fig. 3B) and William plot (Fig. 3C) provided a visual representation of the excellent performance of the best model. As shown in Fig. 3B, the predicted and experimental pIC₅₀ values of modelling compounds exhibited a good linear relationship, indicating good goodness of fit and predictive ability. From the William plot (Fig. 3C), it was found that only compound 30 in the training set had an h value higher than the warning h*, thus being considered a structural outlier. However, its prediction residual was very small and acceptable (<0.12 log units) (see Table S3), also indicating the good predictive reliability of the model. We also found one response outlier (compound 5) with a standardized residual greater than three standard deviation units, thus affecting the model's goodness-of-fit to some extent.

Close

Fig. 3

(A) The weak inter-correlation of nCb- and MDEC-12; (B) Scatter plot of predicted and experimental values of modelling compounds in the best QSAR model; (C) Williams plot of the best QSAR model.

Export to PPT

3.2 Mechanistic interpretation of the modelled descriptors

The selected descriptors in the model equation reflect the role of specific structural and physicochemical properties of the compounds in controlling the response endpoint. These properties will be crucial for the development of potential Nsp14 inhibitors. In line with OECD principle 5 (“a mechanistic interpretation, if possible”) (OECD, 2007), we provided a detailed mechanistic explanation for both descriptors appearing in the best QSAR model, including nCb- and MDEC-12.

Fig. 4A showed the variable importance plot (VIP), from which we found that the most important descriptor in the model (VIP score ≥ 1.0) was MDEC-12 and the second important descriptor (VIP score < 1.0) was nCb-. The loading plot was also used to assess the influence of the two descriptors on the inhibitory activity of Nsp14. The distances between the descriptors and the Y-response or the origin were taken into account, respectively. If one descriptor is close to the Y-response and far from the origin, it is considered more important (Kumar et al., 2020). Thus, it was obvious that from the Fig. 4B, MDEC-12 showed a greater effect on the inhibitory activity against Nsp14, which was consistent with the VIP results.

Close

Fig. 4

(A) Variable importance plot (VIP); (B) loading plot for the two descriptor variables in the best QSAR model.

Export to PPT

Both two descriptors included in the best model showed a positive correlation with Nsp14 inhibitory activity, implying that the higher the nCb- and MDEC-12 values of a compound, the more significant the inhibitory potency on Nsp14. The nCb- descriptor refers to the number of substituted benzene moieties present in the molecule. As shown in Fig. 5, compound 30 had the highest nCb- value (5) in the modelling dataset, indicating that compound 30 has five substituent groups on the benzene ring. Therefore, it exhibited a high pIC₅₀ value (7.523) and was the most active compound in inhibiting Nsp14 in the modelling dataset. On the other hand, the nCb-value of compound 5 (3) was smaller than that of compound 30, implying only three substituent groups on its benzene ring, leading to a relatively decreased pIC₅₀ value (5.301) (Fig. 5). The second descriptor MDEC-12 represents the molecular distance edge between all primary and secondary carbons. As shown in Fig. 5, compound 28 had the largest MDEC-12 value (2.915), thus it showed high inhibitory activity with pIC₅₀ value of 6.943. It is interesting to note that the nCb-value of compound 28 is exactly equal to that of compound 5, both being 3. However, in terms of the MDEC-12 value, compound 28 (2.915) exceeded than compound 5 (1.524), thus, it was observed that compound 28 (pIC₅₀ = 6.943) had a higher pIC₅₀ value than that of compound 5 (pIC₅₀ = 5.301). In the modelling dataset, the pIC₅₀ value (4.85) of compound 1 was the smallest, which was due to both the two descriptors had very small values (nCb-=2; MDEC-12 = 0). The mechanistic interpretation provided a more transparent and understandable relationship between the structural features of a compound and its inhibitory effect on Nsp14. For more information on the descriptors of the compounds, please see Table S3.

Close

Fig. 5

Schematic diagram of the selected descriptors in controlling the inhibitory activity against Nsp14 in the best QSAR model.

Export to PPT

3.3 Initial screening using the best QSAR model

We used the best QSAR model to predict the true external compounds without experimental values, with the aim to look for compounds that might have a significant inhibitory effect on Nsp14 activity. For the 263 untested compounds collected from the PubChem website (https://pubchem.ncbi.nlm.nih.gov/), the best model was used to predict their inhibitory activity against Nsp14. Insubria plot, constructed by the leverage values versus the predicted toxicity, can show whether the true external compounds are within the AD range of the model (Gramatica et al., 2013; Gramatica et al., 2014). As shown in Fig. 6A, the coverage of the model for 263 compounds did not seem to be acceptable, reaching only 47 %, with a significant proportion of compounds being outside the AD range. Therefore, we further screened the true external set and selected 42 structurally most similar compounds from the 263 compounds to as a small true external set (see section 2.6). Obviously, the model showed a significant coverage (81 %) for the 42 compounds (Fig. 6B). At the same time, to further validate the prediction results given by the best model, the prediction quality of each compound was evaluated using the “prediction reliability indicator” tool developed by Roy's team (Roy et al., 2018). For the 263 untested compounds, most were “good” or “moderate” (>90 %), while 26 compounds were “bad” (Table S4). Among the small true external set consisting of 42 untested compounds, 39 compounds were “good” or “moderate” (>92 %) and only 3 compounds were “bad” (Table S5). Therefore, the best model we developed can predict the inhibitory activity of compounds without experimental values to a certain extent. If the structural similarity between the query compounds and the source compounds is better, the predicted results should be more reliable.

Close

Fig. 6

(A) Insubria plot of 263 untested true external compounds; (B) Insubria plot of 42 untested, structurally most similar compounds.

Export to PPT

3.4 Molecular docking analysis

The 42 structurally most similar compounds were docked to the N7-MTase active pocket of Nsp14, of which 17 molecules were successfully docked, and detailed molecular docking analysis was reported in Table S6. Based on the binding energy and predicted inhibitory activity, the top three molecules were evaluated together, including 164775752, 164775726, and 11454312. The predicted pIC₅₀ of compound 164775752 was 6.120, and the binding energy was −29.761 kcal/mol. The predicted pIC₅₀ of 164775726 was 6.432, and the binding energy was –23.199 kcal/mol. Obviously, the comprehensive performance of these two compounds was good. In contrast, for compound 11454312, although its binding energy (-18.219 kcal/mol) was not very satisfactory, its predicted inhibitory activity (pIC₅₀ = 6.767) was the second highest among these 17 molecules, thus it was used as the third target molecule for further study. Although compound 125506109 exhibited the highest predicted activity (pIC₅₀ = 6.9923), it possessed an undesirable binding energy (-8.626 kcal/mol), thus was not selected for further analysis. As shown in Table 2, a detailed summary of the binding energy (CDocker energy), interacting residues, forces involved, and predicted activities of the three compounds was available. At the same time, we also discussed the detailed docking interactions of the selected three compounds with the active site residues of Nsp14.

The best docking conformations of the three compounds 164775752, 164775726, 11454312 with Nsp14 were shown in Fig. 7. The 3D figures (left) highlighted the classical hydrogen bonds (pink) formed between the compounds and amino acid residues. The three compounds also formed non-classical hydrogen bonds and other bonds with amino acid residues, these interactions were shown in the 2D figures (right). The compound with the best binding energy in the screening list was compound 164775752 (Fig. 7A), which interacted with residues ASN306, PRO335, ARG310 through hydrogen bonding, including one conventional hydrogen bond with ARG310 (pink), two carbon hydrogen bonds with ASN306 and PRO335. The compound also formed π-Sulfur bond with CYS309, attractive charge with ARG310, π-π Stacked with PHE426, and π-alkyl hydrophobic bond with LYS423.

Close

Fig. 7

Molecular docking interactions of compounds (A) 164,775,752 (B) 164,775,726 and (C) 11,454,312 with the active site of Nsp14 N7-MTase domain.

Export to PPT

The second-ranked binding energy was compound 164775726 (Fig. 7B). It interacted with residues ARG310 and LYS423 through hydrogen bonding, which included three conventional hydrogen bonds with ARG310 and other one with LYS423 (pink). Compound 164775726 also formed π-π Stacked bond and π-π T-shaped bond with PHE 426, and attractive charge with ARG310.

Compound 11454312 had the highest predicted inhibitory activity (pIC₅₀ = 6.7674) although it exhibited the worse binding energy (Fig. 7C). It interacted with residues ARG310, LYS423 and PHE506 through hydrogen bonding, which included three conventional hydrogen bonds with ARG310 and other one with LYS423 (pink). Furthermore, it formed π-π Stacked bond and alkyl hydrophobic bond with PHE426, π-Sulfur bond with TRP292, two π-alkyl hydrophobic bonds with CYS309. ARG310 was also involved in the formation of attractive charge and π-Cation. As well, LYS423 participated the formation of alkyl hydrophobic bond.

The docking analysis of the three top compounds with N7-MTase active site of Nsp14 showed that they could interact with important amino acid residues such as PHE426 and ARG310. Especially the key amino acid residue ARG310, all the three candidates could interact with it by hydrogen bonding interactions. Therefore, the three compounds may be potential Nsp14 inhibitor candidates.

3.5 ADMET profiling of the top selected screened compounds

The three compounds were further tested for their ADMET properties. As can be seen in Table 3, the ADME properties of all three compounds fell within the acceptable range (O'Brien and Fallah, 2013). Although all three compounds violated one principle (the number of hydrogen receptors were greater than 10) about “Lipinski’s principle of five”, they were acceptable. Regarding toxicity, compound 164775726 was predicted to be potentially hepatotoxic, thus it was excluded. The other two compounds were predicted without significant toxicity, suggesting that they are amenable to being pursued as therapeutic candidates. More details of the ADMET results were provided in the Table S7. To date, no specific drug targeting SARS-CoV-2 Nsp14 on the market, we were unable to find available data to use as references. As a result, compounds 164775752, 11454312 as candidate Nsp14 inhibitors that passed the ADMET predictions need to undergo a series of rigorous clinical trials.