Searching for anthranilic acid-based thumb pocket 2 HCV NS5B polymerase inhibitors through a combination of molecular docking, 3D-QSAR and virtual screening

Abstract

A combination of the following computational methods: (i) molecular docking, (ii) 3-D Quantitative Structure Activity Relationship Comparative Molecular Field Analysis (3D-QSAR CoMFA), (iii) similarity search and (iv) virtual screening using PubChem database was applied to identify new anthranilic acid-based inhibitors of hepatitis C virus (HCV) replication. A number of known inhibitors were initially docked into the “Thumb Pocket 2” allosteric site of the crystal structure of the enzyme HCV RNA-dependent RNA polymerase (NS5B GT1b). Then, the CoMFA fields were generated through a receptor-based alignment of docking poses to build a validated and stable 3D-QSAR CoMFA model. The proposed model can be first utilized to get insight into the molecular features that promote bioactivity, and then within a virtual screening procedure, it can be used to estimate the activity of novel potential bioactive compounds prior to their synthesis and biological tests.

• 3D-QSAR CoMFA
• HCV replication inhibitors
• molecular docking
• PubChem database
• virtual screening

Introduction

Hepatitis C virus (HCV) constitutes a serious worldwide health concern that impacts approximately 180 million people. Approximately, 60–80% of HCV infections evolved in chronic hepatitis usually lead to liver fibrosis, cirrhosis and hepatocellular carcinoma1^,2. Hitherto, the current standard of care (SOC), which combines pegylated alpha interferon and the nucleoside analogue ribavirin, has limited efficacy (approximately 50% of the treated patients achieved sustained virologic response, SVR) and may be associated with a number of side effects3^,4.

In order to address this need, new direct-acting antiviral agents (DAAs) are developed that target viral proteins critical for HCV replication. Recently, the NS3/4 A protease has become a target and two DAAs, boceprevir and telaprevir, in combination with the SOC have been approved as treatment for hepatitis C5^,6. Although these new triple combination regimens improved SVR rates and shortened duration of therapy to HCV genotype 1-infected patients7^,8, significant side effects, such as a low barrier to resistance, appeared. Up to date, several antiviral compounds including NS5A inhibitors, nucleoside/nucleotide analogues and non-nucleoside inhibitors (NNIs) that target the NS5B polymerase, are in clinical development.

HCV, a member of the Flaviviridae family, is a small, single-stranded, positive-sense RNA virus. Its genome encodes a single polyprotein of 3010 amino acids with four structural proteins located in the amino terminus and six non-structural proteins encoded in the remainder. On the basis of the genomic variability in a small region of the gene that encodes the viral NS5B polymerase, the naturally occurring variants are classified into six major genotypes (genotypes 1–6) and several subtypes9.

NS5B is the RNA-dependent RNA-polymerase (RdRp), which is responsible for the complete copy of the viral RNA genome of HCV. Due to its significant action, NS5B has been assigned extensive structural and biochemical characterization10^,11, and several crystallographic structures of polymerase have been structurally elucidated. It has been shown that the protein adopts a morphology similar to that of a right hand with “fingers”, “palm” and “thumb” subdomains, common to several polymerases12. The search for inhibitors of HCV NS5B RdRp has led to investigation of nucleoside and nucleotide inhibitors (NIs) and non-nucleoside inhibitors (NNIs).

NIs are modified nucleosides and nucleotides that can be bound at the active site of polymerase to compete with natural substrates, acting as non-obligate chain terminators. Importantly, due to the conservation of the active site, the efficiency of nucleos(t)ide inhibitors is similar across all HCV genotypes and isolates13. On the other side, NNIs are chemically diverse small molecules that bind to allosteric sites of the enzyme, probably inducing conformation changes in the RdRp, thus hindering the polymerase from functioning effectively and leading to inhibition of the initiation step in RNA replication. Virtual screening (VS) has already proven its value as a useful tool to search for compounds within large chemical libraries, which are most likely to bind to a drug target. This computational tool has been emerged as a reliable, cost-effective and time-saving technique for the drug discovery. VS has already played a significant role in the discovery of several compounds that have either entered clinical trials or even reached the market14–16.

There are two broad categories of VS techniques: (a) ligand-based virtual screening and (b) structure-based virtual screening17–22. For the first category, a model of a receptor can be built, exploiting the common information (pharmacophore or descriptor values) of structurally diverse ligands that bind to the receptor. Then this model can be used to retrieve other potentially active molecular scaffolds. For the second category, the protein target is available, and the candidate ligands are docked into the target and scoring functions are applied to estimate the likelihood that the ligand will bind to the target with high affinity. To date, several studies have been reported in the literature with the objective of identifying novel potent HCV inhibitors using different in silico protocols and workflows23–34.

In this study, a combination of computational methods of (i) molecular docking, (ii) 3-D quantitative structure activity relationship comparative molecular field analysis (3D-QSAR CoMFA), (iii) similarity search and (iv) virtual screening is used to identify new inhibitors of HCV replication. Efforts were carried out to understand the structural characteristics that affect the binding of 53 anthranilic acid-based analogues with HCV NS5B GT1b polymerase serving as receptor35^,36. These compounds have been synthesized and tested as inhibitors of HCV replication, and some of them have been co-crystallized with HCV polymerase, too. In our work, all the compounds were docked into the Thumb pocket 2 (TP-2) of enzyme, and the docking pose of each compound was subsequently used in a receptor-based alignment for the generation of the CoMFA fields. A computational workflow was followed in order to build a 3D-QSAR CoMFA model and is graphically depicted in Scheme 1.

Scheme 1. Computational workflow used in our studies for identification of new compounds – potent HCV inhibitors.

Methods

As aforementioned, a computational workflow has been developed to study potent anthranilic acid analogues as inhibitors of HCV replication in order to predict the inhibitory activity of analogues included in PubChem database. First, a dataset of 53 anthranilic acid-based inhibitors was docked into the TP-2 of the HCV NS5B polymerase for calculating the docking score. Then, the docking pose of each inhibitor was used in a receptor-based alignment for generating the CoMFA fields. A 3D-QSAR CoMFA model was subsequently built to accurately estimate the activity and, finally, within a virtual screening framework anthranilic acid-based analogues included in PubChem database were prioritized for screening.

Dataset

For the development of the 3D-QSAR model, the original dataset consisted of 53 anthranilic acid-based inhibitors (Table 1)35^,36. These analogues were all tested for inhibition of HCV replication in a biochemical assay using a genotype 1b Δ21NS5B protein35^,36. The 3D structures of anthranilic acid-based analogues were built and adjusted using the Maestro 9.3 (New York, NY) molecular builder. All the hydrogen atoms were added, and the molecules were submitted in full structure optimization using the minimization procedure of MacroModel 9.937. For the minimization, water was selected as solvent and a standard molecular mechanics energy function (OPLS_2005 force field38) and the Polak-Ribiere conjugated gradient method (500 iterations with gradient 0.01 kJ/mol Å)3 were used. The minimized structures were subsequently prepared using LigPrep 2.339, where all the structures were checked to be correct and ionization states were generated.

Table 1. Structures of 53 anthranilic acid-based analogues.

*Numbering of compounds reported in literature35^,36.

Preparation of the receptor

There are five co-crystal complex structures of anthranilic acid-based inhibitors available in Protein Data Bank (PDB), with codes 4J08, 4JU3, 4JU4, 4JU7 and 4JJS. Since no significant conformational differences were observed on the crystal structures, the crystal structure with the highest X-ray resolution (PDB ID: 4JU3, 2.00 Å X-ray resolution) was chosen for the docking and was prepared using the Protein Preparation Wizard implementation in Schrödinger suite 201240. The bond orders and disulfide bonds were assigned, all the hydrogen atoms were added and all the water molecules that separated by a distance of more than 5 Å from the allosteric site of TP-2 were deleted. The Epik 2.0 implementation was used to predict ionization and tautomeric states of the ligand het groups40^,41. The hydrogen-bonding network was optimized by reorienting the hydroxyl groups, amide groups of Asn and Gln residues and by selecting appropriate states and orientations of the imidazole ring in His residues. Finally, using the “impref utility” and the OPLS_2005 force field42, the hydrogen atom positions were optimized by keeping all the heavy atoms in place.

Molecular docking

The Glide implemented in Schrödinger was utilized for the molecular docking procedure40. Glide has shown a good reproducibility of co-crystal ligand conformation, which has also been recommended in literature for its accuracy in molecular docking and scoring43. The binding region was defined by a 12 Å ×  12 Å ×  12 Å box centered in the centroid of the crystal ligand, and the option “docking ligands similar in size to the workspace ligand” was chosen. All other parameters were kept as default. The best pose was the output on the basis of the Glide score and the protein–ligand interactions. The extra precision (XP) mode and the XP GScore were used for the docking and scoring, respectively.

CoMFA analysis

In this study, a “protein-based” alignment was used (the conformations of the studied compounds used for the alignment were chosen by docking these compounds in the TP-2 of RdRp). The molecular alignment is the most sensitive factor that has a significant impact on the 3D-QSAR models. Atomic charges for the aligned molecules were calculated using the Gasteiger–Hückel method, which is a combination of two other charge computational methods: the Gasteiger–Marsili44 method to calculate the σ component of the atomic charge and the Hückel45 method to calculate the π component of the atomic charge. The total charge is the sum of the charges calculated by the two methods. The CoMFA fields are generated by creating a grid around the molecule and calculating the steric and electrostatic potentials at each point on the grid using a charged probe atom. The CoMFA calculations were performed using Tripos Advance CoMFA46 module in SYBYL 8.0. The steric and electrostatic field energies were calculated using the Lennard-Jones and the Coulomb potentials, respectively, with a 1/r² distance-dependent dielectric constant in all intersections of regularly spaced (0.2 nm) grid. A sp3 carbon atom with radius of 1.53 Å and charge +1.0 was used as a probe to calculate the steric and electrostatic energies between the probe and the molecules using the standard Tripos force field47. The truncation for both the steric and electrostatic energies was set to 30 kcal mol⁻¹. This indicates that any steric or electrostatic field value that exceeds this value will be replaced with 30 kcal mol⁻¹, thus makes a plateau of the fields close to the center of any atom. An optimal number of components (ONC) was designated such the cross-validated r² was maximal, and the standard error of prediction was minimal.

Partial least square analysis

The partial least square (PLS) approach applied in SYBYL 8.048^,49 was used to derive the 3D-QSAR model and linearly correlate the CoMFA fields to the pIC₅₀ values. The CoMFA fields were used as independent variables and the pIC₅₀ values as dependent variables. PLS was performed in two steps. In the first step, the cross-validation analysis was performed using Leave-One-Out (LOO) method, which was applied to determine the value of the cross-validated r² (q²), the cross-validated standard error of predictions, S_PRESS and the ONC. In the LOO method, one compound is removed from the dataset, and the model is derived involving the rest of the molecules. Then using this model, the activity of the omitted molecule can be predicted50. The model with the highest q² and the lowest Standard Error Prediction (SEP) was driven for further analysis. The cross-validated r² is calculated by the following equation: (1) where Y_obs and Y_pred indicate observed and predicted activity values, respectively, and indicates mean activity value. A model is considered acceptable when the value of q² exceeds 0.551.

In the second step, using the ONC obtained from cross-validation analysis, a non-cross-validated method was applied and the non-cross-validated R² (conventional) is calculated according to the following equation: (2) In addition, the statistical significance of the model was described by the standard error of estimate (SEE) and the probability value (F-value). To speed up the analysis and reduce noise, a minimum column filter value σ (2.00 kcal mol⁻¹) was used. The application of non-cross-validated method led to the build of CoMFA model.

3D-QSAR model validation

The predictive capability of the 3D-QSAR CoMFA model was evaluated with an external test set using the molecules that have not been taken into account during the process of developing the model52. The predictive correlation () is calculated using the following equation: (3) where Y_{pred (test)} and Y_test indicate predicted and observed activity values of the test set, respectively, and indicates mean of observed activity values of the training set. For a predictive QSAR model, the value of should be more than 0.5.

Furthermore, external predictive power of the developed 3D-QSAR model using test set was examined by considering . The formula is given below53: (4) where r² and are squared correlation coefficients between the observed and predicted activities of the test set with and without intercept, respectively. is a new, rigorous and powerful statistical index reported by Roy et al.54. In their study, they demonstrated that even when the value of q² is high, this is not a sufficient condition to indicate the predictive capability and power of the model. For a significant external model validation, the value of should be greater than 0.5. The index can be calculated using a web application accessible from http://aptsoftware.co.in/rmsquare/. The concept of can also be applied for the training set.

To assess the predictive power of produced 3D-QSAR model, Enalos Model Acceptability Criteria KNIME node was used55 that includes all proposed tests by Tropsha. Especially, for continuous QSAR, criteria followed in developing activity/property predictors are as follows: (i) correlation coefficient R between the predicted and observed activities; (ii) coefficients of determination (predicted versus observed activities and observed versus predicted activities for regressions through the origin); (iii) slopes k and k' of regression lines through the origin. The QSAR model is acceptable if the following conditions are satisfied: (i) q² > 0.5; (ii) R² > 0.6; (iii) and 0.85 ≤ k ≤ 1.15 or  < 0.1 and 0.85 ≤ k′ ≤ 1.15; and (iv) 52^,55.

Domain of applicability

When proposing a validated model, it is very important to simultaneously define its limits so that a well-defined applicability domain could indicate those predictions that can be considered reliable. When the model is used to screen new compounds, it is important that these structures fall out the domain of applicability of the model and are filtered out as the model cannot generate reliable predictions for these structures. Domain of applicability can be defined using similarity measurements based on the Euclidean distances among all training and test compounds. The distance of a test compound to its nearest neighbor in the training set is compared to a predefined threshold (APD), and the prediction is considered unreliable when the distance is higher than that. APD was calculated based on the following formula: (5)

Calculation of <d> and σ was performed as follows: first, the average of Euclidean distances between all pairs of training compounds was calculated. Next, the set of distances that were lower than the average was formulated. <d> and σ were finally calculated as the average and standard deviation of all distances are included in this set. Z is an empirical cutoff value and for this work, it was chosen equal to 0.556–59. Enalos Domain – similarity node that executes the aforementioned procedure is included in our workflow and was used to assess domain of applicability of the proposed model56–59. For the calculation of the domain of applicability, we have used Enalos Domain KNIME nodes. These nodes have been developed by Novamechanics Ltd. (Nicosia, Cyprus) and are publicly available through the KNIME Community and the company’s website60.

Virtual screening

PubChem database has been recently emerged as a very important source of compounds that could be prioritized for screening for a wide range of targets. We have decided to use PubChem database within our in silico framework to virtual screen compounds that could act as HCV inhibitors. For this purpose, we have developed a KNIME workflow to assess the compounds that will be prioritized for virtual screening. First, our purpose was to mine all anthranilic acid-based analogues within PubChem database. To achieve that, we have used Substructure Matcher node included in Indigo cheminformatics toolkit. All these compounds were compared to the most active compound of our original dataset in terms of similarity given by Indigo Fingerprints. It is known that similarity searching refers to the calculation of the similarity coefficients for a given molecule and each compound included within a database. Among the different proposed methodologies, similarity searching using fingerprints is one of the most widely applied methods for ligand-based screening. Fingerprints are representations of the structural features of a molecule given as a series of binary digits (bits) that account for the presence or absence of particular substructures in the molecule. For this study, Fingerprint Similarity node using Tanimoto metric for similarity computation was used to afford the most similar compounds that could be prioritized for screening.

All the retrieved anthranilic acid-based analogues were virtually screened with the ultimate goal of prioritizing the most promising new synthetic targets. The proposed 3D-QSAR model was used to identify the most active compounds that are proposed for further investigation.

Results and discussion

Molecular docking of analogues into the TP-2 of NS5B HCV polymerase

The TP-2 is a narrow hydrophobic cavity near the base of the thumb domain, approximately 35 Å apart from the active site. Sequence analysis showed that the primary structure of this binding site was relatively well conserved among the known HCV genotypes61. The binding site traverses the thumb from front to back of the RdRp “right hand” but it is relatively shallow and narrow (30 Å ×  10 Å ×  10 Å). It is bounded by three of the α-helices in the thumb subdomain separated by the link between two of these from the fingertip α-helix/site I region that is 15 Å away and is proximate to the C-terminus of the thumb. The mechanism of action for the inhibitors of TP-2 is unknown. It has been suggested that inhibition is due to blocking association with other proteins that are important for replication61. The TP-2 consists of a hydrophobic pocket with residues M423, W528, L419, Y477 and R422, while hydrogen bonds are favorable. The observed key interaction between the anthranilic acid-based inhibitors and the enzyme is two hydrogen bonds between the carboxylate oxygens and the backbone amide nitrogens of S476 and Y477.

The molecular docking was first applied to anthranilic acid-based and active inhibitor 8 (ID = 23, Table 1), which it has already co-crystallized with the HCV NS5B polymerase and described by Stammers et al.35. In Figure 1, the superposition of crystal structure conformation and that obtained in the thumb pocket of compound 8 is depicted. Thus, it is obvious that the algorithm Glide used for molecular docking is able to reproduce crystallographic experimental data for structural similar compounds. Afterward, the 53 anthranilic acid-based molecules were docked into the TP-2 allosteric site of HCV polymerase. The docking results expressed in XP GScore are listed in Table 2. The lower the value of XP GScore, the more inhibition ability is expected for the compound. A good correlation was found between the docking scores and the IC₅₀ values.

Figure 1. Superposition of two structures; crystal structure is depicted with green and docked compound 8 with orange (RMSD ∼ 0.2).

Table 2. The IC₅₀ values, the XP GScore, the pIC₅₀, the predicted pIC₅₀ and the residuals for the anthranilic acid-based analogues.

ID	IC50 (μM)	XP GScore	pIC50	Predicted pIC50	Residuals
1	26	−7.22	4.585	4.621	−0.036
2*	9	−7.949	5.046	4.698	0.348
3	11	−6.724	4.959	5.037	−0.078
4	43	−7.462	4.367	4.683	−0.316
5	8.4	−6.395	5.076	5.014	0.062
6*	4.5	−7.772	5.347	4.913	0.434
7	6.1	−6.357	5.215	5.121	0.094
8	6.8	−7.349	5.167	5.025	0.142
9	4.6	−6.139	5.337	5.081	0.256
10	3.7	−6.739	5.432	5.523	−0.091
11	6.7	−6.525	5.174	5.21	−0.036
12	0.55	−6.138	6.26	6.303	−0.043
13	0.22	−7.369	6.658	6.489	0.169
14*	0.34	−7.618	6.469	6.066	0.403
15	0.61	−8.087	6.215	6.317	−0.102
16*	0.16	−8.86	6.796	6.276	0.52
17	0.31	−7.845	6.509	6.419	0.09
18	0.78	−8.237	6.108	6.101	0.007
19	0.47	−7.535	6.328	6.535	−0.207
20*	0.53	−7.632	6.276	6.511	−0.235
21	0.29	−9.146	6.538	6.513	0.025
22	0.33	−6.961	6.481	6.5	−0.019
23	0.085	−12.606	7.071	7.072	−0.001
24	0.27	−8.971	6.569	6.58	−0.011
25*	0.37	−7.468	6.432	6.339	0.093
26	0.21	−7.996	6.678	6.73	−0.052
27	1.9	−6.217	5.721	5.633	0.088
28	1.1	−7.41	5.959	5.98	−0.021
29*	1.5	−7.65	5.824	5.881	−0.057
30	2.9	−6.969	5.538	5.601	−0.063
31	1.4	−7.559	5.854	5.805	0.049
32	2.2	−6.79	5.658	5.657	0.001
33	3	−8.554	5.523	5.521	0.002
34*	1.64	−7.936	5.785	6.102	−0.317
35	0.67	−8.279	6.174	6.093	0.081
36	1.3	−8.036	5.886	6.06	−0.174
37	2	−7.841	5.699	5.622	0.077
38	1.4	−8.148	5.854	5.938	−0.084
39*	0.74	−9.008	6.131	6.257	−0.126
40	4.3	−6.555	5.367	5.553	−0.186
41	1.3	−8.811	5.886	5.776	0.11
42	0.46	−9.069	6.337	6.343	−0.006
43*	>40	−8.42	–	5.607	–
44	0.19	−8.898	6.721	7.005	−0.284
45	1.3	−7.732	5.886	5.882	0.004
46	0.47	−7.496	6.328	6.224	0.104
47	0.74	−8.489	6.131	6.137	−0.006
48*	0.45	−8.602	6.347	6.476	−0.129
49	0.33	−8.219	6.481	6.297	0.184
50	0.3	−9.087	6.523	6.58	−0.057
51	0.2	−8.464	6.699	6.485	0.214
52	0.08	−9.211	7.097	6.994	0.103
53	0.19	−9.026	6.721	6.71	0.011

*Test set.

The molecular docking results of compound 8 of this series are represented in Figure 2. As aforementioned, the pose that adopted is the same with the co-crystal inhibitor. The right hand side (RHS) of anthranilic acid scaffold (4-methylphenylsulfonamide group) is placed into the deep lipophilic pocket of enzyme and a new hydrogen bonding interaction between the left hand side (LHS) of anthranilic scaffold, and the side chain of R490 is added to the three hydrogen bonds with S476, Y477 and L497 appeared in the X-ray crystallography. No significant differences to the positioning and orientation of the anthranilic acid of other analogues were observed. All anthranilic acid-based derivatives, except from compounds 10a and 18, adopt a binding mode to enable bidentate hydrogen bonding between the analogue carboxylic acid and the backbone N-H’s of S476 and Y477. This hydrogen bonding pattern is similar to a reported inhibitor35.

Figure 2. The representative docking results of active analogue 8 with the highest docking score −12.606. The hydrogen bonding interactions are depicted with yellow dotted lines.

Dataset alignment

The present CoMFA model was generated using a “protein-based” alignment of the training set consisting of 42 anthranilic acid-based inhibitors. The results of the alignment are depicted in Figure 3. The active analogue 8 was used as a template molecule for the alignment using 12 atoms (Figure 4) that are common in all analogues. A rigid-body atom-by-atom superimposition of one molecule onto another was performed using the utility “Align database” available in SYBYL 8.0.

Figure 3. Database alignment superposition of training set used for 3D-QSAR analysis based on docked conformation of 8.

Figure 4. Common substructure of the 53 compounds subjected to the alignment.

3D-QSAR CoMFA model

Inhibitory concentrations (IC₅₀) of the 53 anthranilic acid-based inhibitors were converted into corresponding pIC₅₀ (−log IC₅₀) and used as dependent variables in CoMFA analysis. The data set was classified in a way that it could be obtained a representative training set in terms of molecular structure and diverse bioactivity62–77. The distribution of the activity values for the test set follows the distribution of the activity values for the training set. According to Golbraikh and Tropsha78, this approach is correct since representative points of the test set must be close to those of training set and vice versa. The separation of the dataset into training and test set was performed in a commonly used ratio of 4:1. The training set consists of 42 molecules, and the test set consists of 11 molecules (Table 2). The Partial Least Squares (PLS) analysis was used for deriving the relationships among the CoMFA fields and the inhibitory activity (IC₅₀ values)49^,50. The 3D-QSAR CoMFA model was derived using the standard implementation in SYBYL 8.0 molecular modeling package46.

CoMFA analysis

CoMFA models are generated using the combination of steric and electrostatic fields. An ONC of five corresponds to the highest cross validated q² of 0.708 and the lowest S_PRESS of 0.376. A non-cross validated of 0.965, a SEE of 0.130 and an F value of 199.276 with zero probability of show a good statistical correlation between the predicted and the experimental pIC₅₀ values for the non-cross validated form of the CoMFA model. In Figure 5, the relationship between the predicted and the observed pIC₅₀ values is represented. The CoMFA analysis shows that the relative field contributions are 70.4 and 29.6% for steric and electrostatic field, respectively.

Figure 5. Plot of the predicted pIC₅₀ values versus observed ones for the CoMFA model.

Validation of the 3D-QSAR model – domain of applicability

External validation79^,80 was also performed to further assess the stability and the predictive ability of the CoMFA model. This validation was performed using 11 molecules not included in the construction of the model. The predicted pIC₅₀ values are in correlation with the observed ones within the tolerable error range (Figure 5). The external has a value of 0.820 upholding good correlation between predicted and observed pIC₅₀ values. The proposed model has values equal to 0.965 and 0.525 for the training and test set, respectively (Table 3). Moreover, the Enalos Model Acceptability Criteria KNIME node has been applied to the data. The model passed Tropsha’s recommended tests for predictive ability and the results are depicted in Figure 6. These results suggest that this alignment can effectively take into consideration the ligand-receptor interactions, the CoMFA model is thus reliable and could be used in the design of new inhibitors of HCV replication within this structural motif of molecules.

Figure 6. Enalos Model Acceptability Criteria KNIME node screenshot.

Table 3. Statistical parameters of the CoMFA model.

CoMFA model
ONC	5
q^2	0.708
SPRESS	0.376
	0.965
SEE	0.130
F (n1 = 5, n2 = 36)	199.276
Prob. of r^2 = 0 (n1 = 5, n2 = 36)	0.000
Field contribution steric:electrostatic	0.704:0.296
	0.820
(training set)	0.965
(test set)	0.525

It is important that the limitations of the model are also described via the domain of applicability. This gives an important indication as the user can freely and creatively design novel molecules but will be warned for the reliability of the prediction when the structural characteristics cannot be tolerated by the model. After model validation, the domain of applicability of our model was also defined to ascertain that a given prediction can be considered reliable56–59. Enalos Domain – Similarity node that executes the aforementioned procedure is included in our workflow and was used to assess domain of applicability of the proposed model56–59. The applicability domain limit value was defined equal to 19.537 based on the equation provided in Methods section. All compounds in the test set had values in the range of the domain of applicability. The predictions for all compounds that fell inside the domain of applicability of the model can be considered reliable.

CoMFA contour maps

The results of a CoMFA analysis were graphically represented by field contour maps. The StDev*Coeff option was selected to control the type of CoMFA fields to be viewed. The default option “contribution” was selected to control the interpretation of the CoMFA fields. In these maps, the areas around the compound that interacts favorably or unfavorably with the receptor are visualized. The contour maps are used to identify the structural features relevant to the biological activity in this series of analogues and their obtained information can be utilized for further design of anthranilic acid-based inhibitors of HCV replication.

The CoMFA fields around the active compound 8 are displayed in Figure 7. The steric interactions are represented by green and yellow contour maps (Figure 7a); where bulky groups near the green regions increase the inhibitory activity, they cause the opposite effect near the yellow regions. The electrostatic interactions are displayed by red and blue contour maps (Figure 7b); where electronegative groups near the red regions increase the inhibitory activity and they decrease it near the blue regions.

Figure 7. CoMFA StDev*Coeff contour maps around the active compound 8: (a) solid for steric field (green: bulky groups are favored; yellow: bulky groups are disfavored); (b) mesh for electrostatic field (red: electronegative groups are favored; blue: electropositive groups are favored).

As it can be observed, the yellow regions are more than green ones for steric field in the contour map of the CoMFA model (Figure 7a). The large green contour map and small one around the LHS group of the anthranilic acid scaffold indicate that bulky groups at these positions increase the activity. For example, compounds 4, 5, 9, 10 and 40–47 bearing bulky groups, such as –CF₃ on o- or m-position of phenoxy group showed high activities. Similarly, the small green contour map above the atom of nitrogen justifies the high activities of compounds 24, 36, 38 and 40–47, as N-iPr replaces the N-H group. On the contrary, there is a yellow contour map around the RHS group on the anthranilic acid scaffold in the contour map for steric field that indicates low activities of compounds bearing bulky groups on m-position of benzyl ring, such as compounds 5a–c. The second big yellow contour map indicates that the presence of a linker, such as the oxygen linker, between the two rings increases the activity that does not happen with inactive compounds 9a and 10a. Finally, the other three yellow contour maps on the RHS of anthranilic acid scaffold show that bulky groups among the nitrogen atom and the ring are disfavored, while a sulfonamide or carboxamide groups are suitable for this position.

In the contour map for electrostatic field in Figure 7(b), there is a huge blue region around the binding site indicating that the introduction of electropositive groups around this position would increase the inhibitory activity. Nevertheless, the big red contour map and the small one around the LHS group of the anthranilic acid scaffold indicate that compounds with electronegative groups in these positions possess higher activity. Considering the positive potential of R490 in the proximity of LHS group, it will be suitable for negative groups, such as carboxyl groups, to be placed in this position (compounds 7 and 8).

Virtual screening

As described in the Materials and Methods section, the produced QSAR model can be used as a useful tool for screening existing databases or virtual libraries in an effort to prioritize potent compounds for experimental evaluation. In this context, a virtual screening of anthranilic acid-based analogues that are included in PubChem database was carried out using a KNIME workflow. First, 160112 compounds containing the anthranilic acid scaffold were retrieved from the PubChem database. These compounds were compared to the active compound within our original database, compound 8 (ID = 23, IC₅₀ = 0.085 μΜ), using fingerprints based on the Tanimoto similarity metric as described in Materials and Methods section. The 1279 compounds with a Tanimoto similarity metric of over 0.80 were selected for prioritization based on the prediction of the proposed CoMFA model that was used to assess the potency of these compounds.

Our approach that is based on the scaffold and fingerprint similarity search within PubChem database resulted in a narrower chemical space that increases the chance of success since our 3D-QSAR was built based on the specific scaffold.

The pIC₅₀ values of the 18 most potent compounds (1vs–18vs) predicted with the established CoMFA model are listed in Table 4. All the virtual screening compounds were within the domain of the applicability of the model, and therefore their activity estimations can be accepted with confidence. The predicted compounds were searched for their activity using the PubChem and ChEMBL databases, but no biological data has been found. Our searching was continued in SureChem patent chemistry database that contains a collection of over 15 million chemical structures of United States (US), European Patent (EP) and World Intellectual Property Organization (WIPO) full-text patents since 1976. The information retrieved is summarized in Table 5.

Table 4. Structures and pIC₅₀ values of VS compounds.

Table 5. Documents of VS compounds based on SureChem patent chemistry database.

VS compound	Publication date	Assignee/applicant	Description
1vs	24 December 2008	Tibotec Pharmaceuticals Ltd.	A homogeneous time resolved fluorescence based test system
	11 October 2007
2vs–14vs	30 September 2008	Tibotec Pharmaceuticals Ltd.	Small molecule entry inhibitors
15vs–18vs	No data has found

It is notable that compounds 2vs–14vs have been introduced as small molecule entry inhibitors for HIV. In the patent no.: US7429677B281, it is reported that these compounds are inhibitors of the entry process of the HIV into the host cell. More specifically, they are drugs designed to block HIV from entering the human cell by interfering with various phases of attachment and fusion between HIV and the cell. As the co-infection of HCV and HIV has become a worldwide public health challenge, it is very important that our proposed methodology managed to retrieve and to conclude in these small molecules that have been already proven to act as HIV inhibitors as these small molecules could be potential drug candidates as dual inhibitors of HIV and HCV infection.

Conclusions

In this work, we have successfully developed a computational virtual screening workflow combining Molecular Docking, 3D-CoMFA QSAR and similarity search that was applied to PubChem database with the aim to identify anthranilic acid based analogues that could act as inhibitors of HCV replication. In the proposed approach, we have first developed a robust, validated and predictive 3D-QSAR model that was subsequently used to virtually screen compounds from PubChem database. We have narrowed the chemical space search by focusing only on compounds that contain the anthranilic acid scaffold. Similarity search was then used to retrieve the compounds that are similar to known active inhibitors. In this way, we have managed to identify the most promising compounds from a pool of new analogues and to prioritize a list of compounds for screening. The approach revealed several promising chemistry driven compounds with potential high activity. The workflow can also be used to screen other databases or virtual combinations to identify derivatives with desired activity.

Declaration of interest

Authors report no declarations of interest.

This work is supported by funding under the Seven Research Framework Program of the European Union. Project SYSPATHO (ex. PATHOSYS – HEALTH-F5-2010-260429).