Docking and 3D-QSAR Studies on Some HCV NS5b Inhibitors

A theoretical study has been carried out to interpret and support experimental findings regarding inhibition mechanism of HCV NS5b. Twenty-five HCV NS5b inhibitors were docked by QM-Polarized Ligand Docking (QPLD) technique. The comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) methods were used to derive 3D-QSAR models for the selected inhibitors. The CoMFA and CoMSIA models show good cross-validated (Q) and non-cross-validated (R) coefficients for the suggested inhibitors of 0.43, 0.98 and 0.65, 0.99, respectively. The inhibition mechanism was explored and validated. Details of the interactions between the inhibitors and HCV NS5b are given in terms of steric, electrostatic, hydrophobic, hydrogen bonding fields. Enhancing potency via substitutions at positions, which were explored based on these parameters. A good correlation was found between 3D-QSAR and docking results.


Introduction
Hepatitis C Virus (HCV) is a genus of the Flaviviridae family with six major genotypes. It includes a large number of subtypes within each genotype [1]. Worldwide distribution of HCV genotypes includes genotype 1 (Japan, Europe and North America), genotype 2 (Japan and North America), genotype 3 (Indian subcontinent), genotype 4 (North Africa and the Middle East), genotype 5 (South Africa) and genotype 6 (South East Asia) [2,3].
HCV epidemiology is approximately 3% prevalent worldwide [4,5], one-fifth of the globally HCV carriers live in the Middle East [2]. HCV infection was estimated to be about 20% of Egyptians with a newly half million infections per year [6][7][8]. Egypt has the highest globally incidence of HCV infection [9][10][11][12][13] and it varies from 9 to 50% in certain rural areas due to the specific modes of infection [5]. While the mechanism of infection is fully acquainted, the mechanism of viral entry and replication are not completely understood [8]. HCV consists of different enzymes, though there are many clues that HCV polymerase is the maestro enzyme in viral replication process [14,15].
Similar to any DNA polymerases, HCV NS5b adopts the architecture of a right hand with "thumb," "palm," and "fingers" domains. "The palm domain catalyzes the phosphoryl transfer reaction, whereas the fingers domain participates in interactions with the incoming nucleoside triphosphate as well as the template base to which it is paired [14]. However, the thumb domain plays a role in positioning the duplex DNA and in processivity and translocation [14]. HCV NS5b is an interesting target for antiviral therapy with limited side effects. It was the subject of extensive trails to design both nucleoside and non-nucleoside inhibitors [15,16]. Different studies discussed various potential sites of HCV NS5b: Thumb pocket-I (Pro495, Pro496, Val499), Thumb pocket-II (Leu419, Met423), Palm pocket-I (Asn411, Met414, Tyr448) and the active site (Ser282) [17][18][19][20][21][22][23][24].
Molecular docking represents the basics of rational drug design [25], since it predicts the appropriate structure conformation for the potential target. The first docking study was done by Levinthal et al. [26] to predict the conformations of hemoglobin fibers.
Quantitative Structure Activity Relationship (QSAR) was found to be a good tool for reproducing activity of various inhibitors toward particular receptors [27]. The threedimensional QSAR (3D-QSAR) is a good representative of inhibitor-receptor interactions with their bioactive conformations [28,29].
The current work aims at studying some HCV NS5b inhibitors and their interaction with four potential sites (Thumb pocket-I, -II, Palm pocket-I and the active site) using various docking protocols and constructing 3D-QSAR models, such compounds were selected for the present work because of their potency, diversity and selectivity besides the high resolution of the experimental complex structure (2HAI). This investigation helps in understanding of inhibition efficiency of the selected compounds and tailoring of new inhibitors.

Methodology
A set of twenty five inhibitors of diverse activities was selected from a previous experiment [18]. These inhibitors are displayed in Figure 1. They were selected based on their structure diversity, hydrogen bonding ability, hydrophobic moieties, aromatic rings, substitution patterns and their potency toward HCV NS5b allosteric pockets. All inhibitors were sketched by Epike v. 2 module [30]. Then, stereoisomerism and tautomerization options were selected with OPLS-2005 force field. A set of 88 conformers were produced with Ligprep module embedded in Schrodinger suite. These ligands were docked by QPLD method.
QPLD includes Glide and QSite packages. The QPLD algorithm launches with Glide docking job which produces several protein-ligand complexes. After that, QSite performs a single-point energy calculation on each complex treating the ligand with ab initio methods and deriving partial atomic charges using electrostatic potential fitting as QM/MM environment. Glide then re-docks the ligand using each of the ligand charge sets calculated by QSite. Finally QPLD algorithm gives the most energetically favorable pose. The fully automated algorithm is calibrated to provide useful settings according to the user's need [31,32]. The mixed QM/MM approach in QPLD is more accurate especially in the case of metalloproteins [31]. The accuracy in determining binding modes from Glide is about 82% [33] which makes it an appropriate tool for the current docking study. Glide standard error was found to be ~2 kcal/mol.
Before docking, all missing atoms were adjusted to match the original protein structure. Protein preparation wizard was used to assign all bond orders, hydrogen atoms, formal charges on the metallic centers. The neighboring atoms were also assigned and all water molecules were deleted. Prediction of ionization and tautomeric states of the metal groups at pH 7 was performed. A minimization with 0.1 Å RMSD and OPLS-2005 force field was applied. Finally, the whole complex, 2HAI, was split into protein and inhibitor. All inhibitors were docked using Glide XP score in QPLD [34]. Docking validation can be done by calculating coefficient correlation between IC 50 [35,36] or -log IC 50 (pIC 50 ) [37][38][39] and docking scores.
All conformations were analyzed using LPC package [40] to determine the best interaction. The best docked conformations for all inhibitors docked by QPLD were also used for Autodock (Genetic algorithm) and Glide XP calculations. Unfortunately, the scores derived from both methods showed no good correlation with the experiment, see supporting material. Molegro molecular view [33] and Chimera [41] were used to extract the docked inhibitors from their receptors and for their graphic representation, respectively.
The 3D-QSAR models were performed for the most active conformations extracted from docking. These 3D-QSAR models were derived from comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA) methodologies with PLS algorithm [42]. The 3D-QSAR calculations were performed using Sybyl v. x2 which is based on the Tripos force field with a distance-dependent dielectric, the Powell conjugate gradient algorithm and Gasteiger-Hückel charges [43].
CoMFA represents a 3D cubic grid with 2 Å spacing and 4 Å automatic extension away from the investigated molecules in all three axes(X, Y, Z directions). Lennard-Jones and Columbic potentials were calculated to determine the steric and electrostatic interaction fields at each grid intersection. An sp3 carbon atom with a radius of 1.52 Å bearing +1 charge was used as a probe atom to calculate the CoMFA fields [44]. CoMFA contour maps were plotted to analyze steric and electrostatic features to shed light into the strength of binding between the inhibitors and receptor as determined from docking.
CoMSIA calculates five similarity descriptors namely steric, electrostatic, hydrophobic, hydrogen bond donor and acceptor [45]. Gaussian type function was employed to calculate similarity indices at all lattice points. Smoothening function applied by a default value of 0.3. The probe atom with charge +1, hydrophobicity +1, and hydrogen bond donor and acceptor property of +1 was used. The steric indices were related to third power of the radii, partial atomic charges were used to derive the electrostatic fields. Atom based parameters computed hydrophobic descriptors and rule based method was used to get the hydrogen bond donor and acceptor fields [46].

Docking
It was reported that the inhibition of the allosteric sites reduces the catalytic activity of the HCV NS5b [15]. Therefore, there is an unavoidable conflict in rational design of drugs targeting these sites [53]. Li et al. [18] and Love et al. [53] designed a lead for these sites. As shown in Figure 2, the blue part of 1 tends to bind with Thumb pocket-I, while the purple area binds to Thumb pocket-II. The red region binds to Palm pocket-I. For rationalizing the experimental findings, QPLD was applied to elucidate the interaction of dihydropyrones derivatives with Thumb pocket-I, -II, Palm pocket-I and the active site.
Docking validation was done by redocking of inhibitor 19 with RMSD of 0.9 Å when compared to its crystal structure, 2HAI. This assumes a good docking [54], see Figure 3. QPLD gives a good correlation between the calculated binding energies and IC 50 as emerges from R 2 values of 0.68→0.80 and 0.82→0.90 with Thumb pocket-II and Palm pocket-I, respectively. These values are better than that recorded for Thumb pocket-I, Table 1. This indicates that most compounds tend to bind strongly with Thumb pocket-II, Palm pocket-I rather than with Thumb pocket-I. The investigated inhibitors do not show any activity toward the active site as listed in the supporting information. The inhibitor 1 was found to interact selectively with Thumb pocket-II which does not exist in the thumb domain of other polymerases (cellular or viral) [55]. The present study agrees with this view where 1 shows a considerable activity toward all allosteric sites rendering it as a lead compound (IC 50 =0.93 µM), Figure 3. Docking study in the thumb domain, Thumb pocket-I, indicates that the phenol group in compounds 1 and 2 are directed toward the hydrophobic regions (Val485, Leu489 and Val494). Due to the proximity of Thumb pocket-I and -II, they interact similarly with compounds 3-20 as these interactions will be discussed later. However, the rest of compounds do not show a remarkable activity for Thumb pocket-I, besides compound 25 was excluded from the current docking work because of its tiny IC 50 which can't compared with the others, Table 1. For Thumb pocket-II, docking results match the experimental study of Love et al. [55]. As displayed in Figure  3, the inhibitor 1 when compared to crystal structure, 1OS5, the cyclopentyl group lies in hydrophobic residues (Leu419, Arg422 and Met423), the phenol group is surrounded by Leu497, Ile482. Besides, there are two hydrogen bonds formed by the diketo groups of dihydropyran-2-one with Arg501 and Ser476. Compounds 1-16 react with amino acids residues in a similar manner except 2, 3, 6, and 12 which show little differences, see Figure 3.  . exhibit Crystal structure of (2HAI) compared to docked structure of inhibitor 19, crystal structure (1OS5) compared to docked structure of inhibitor 1. Thumb pocket-I, Palm pocket-I poses for inhibitor 1, 1-16 except for 2, 3, 6, 12 and 2, 3, 6, 12 also displayed next to it respectively. Amino acid residues appear as thin sticks while ligand atoms are represented as bold sticks. The hydrophobic parts appear in green while the hydrophilic moieties in orange and gray colors. Hydrophilic residues have a red color, while hydrophobic residues have blue color. Atoms of residues are colored according to the hydropathy index proposed by Kyle and Doolittle in 1982 (see http://en.wikipedia.org/wiki/Hydropathy_index for details). The blue dashed line represents the hydrogen bond.
The phenol group of compound 2 is oriented toward Arg422 and surrounded by Tyr477 and Leu419, parallel stacked with Trp528, hydrogen bonded with Ser476 and Arg422 through O2 and O4, respectively. This is confirmed by unfavorable hydrogen bond acceptors (deep red polygon for O2) steric and hydrophobic (yellow and gray for the phenol group) unfavored region, Figure 4. However, compound 2 doesn't fit well with this pocket because of its small size. Therefore, potency of this compound is decreased by 100 times compared to 1. For compound 3, the anisole group is located close to the far side of Thumb pocket-II, Val485, Leu489, while there are two hydrogen bonds formed through O3 and O2 atoms of the diketo groups of dihydropyran-2-one with Ser476 and Arg501. This inversion of the structure might explain the double potency of 3 compared to 2 (51 vs. 100 µM). For compound 6, three hydrogen bonds are formed between Ser476, Arg501, Lys533 and O2, O3 and O4, respectively, which is explained by unfavorable hydrogen bond acceptors (deep red polygon for O2, O3) and disfavored hydrogen bond acceptor (magenta contour for O4), see Figure 4. This slightly increases the potency (IC 50 =52 µM) compared to 3. The diketo groups of dihydropyran-2-one in 12 form hydrogen bonds with Ser476 and Arg501 through O2 and O3. Moreover, the toluene moiety interacts with the hydrophobic region (favorable/orange and grey/unfavorable) of Leu419, Met423, Ile482, Val485 and Leu489. This decreases the relevant IC 50 to 23 µM.
Cirrhotic patients or patients with decompensated liver are difficult group to cure, those patients with unmet medical needs will deteriorated through 4 years when they treated with the current approved drug [56].
Due to financial issues, we cannot implement any experimental work on the compounds under study, but using the current docking study as a test set to screen the drug database for isosteres [57], which can be used as potential HCVNS5b drugs for cirrhotic patients [56] in another future work. Increasing the size of hydrophobic groups at the para position leads to higher potency (for compounds 6, 8→10 where IC 50 decreases from 52, 9.7 to 1.7 µM; green, orange polygon sterically favored hydrophobic interaction). Further increase in potency is achieved with addition of a fluorine atom at the meta position of the phenyl ring (for compounds 17-19, IC 50 decreases from 0.89 to 0.53 µM; orange favored hydrophobic interaction exist in both meta positions).
The bioactive conformations for Thumb pocket-II were aligned and pharmacophore features were mapped by PharmaGist [58]. Figure 5 displays inhibitor 1 docked in Thumb pocket-I, -II and Palm pocket-I. Hydrophobicity (green), hydrophilicity (orange), aromaticity (red) and nonpolar hydrogen bonds (white) of 1 are shown in Figure 3. The pharmacophore models of bioactive conformations for Thumb pocket-I, -II suggest one aromatic center (phenol group as given in pink color in Figure 2), three hydrophobic centers (ethyl group (in black and the chiral carbon of pentane ring (blue)) and two acceptors (two carbonyl groups of dihydropyran-2-one ring (blue)). The aromatic center interacts with the side chains of Tyr 477 for compounds 1 and 2. However, there is no such interaction in other inhibitors. The hydrophobic centers interact with Leu419, Arg422 and Tyr477. Most inhibitors form hydrogen bonds with Ser476, Arg501 residues, except 2, 22 which form hydrogen bonds with other residues. These results coincide with docking data with little differences.
For Palm pocket-I, the pharmacophore model of the bioactive conformations shows three acceptors (two carbonyl groups and oxygen atom of dihydropyran-2-one ring(blue)) and three hydrophobic centers (two ethyl groups (black) and the chiral carbon of pentane ring (blue)). The three acceptors form hydrogen bonds with Leu446, GLY449, while third residue varies from one inhibitor to another. The hydrophobic centers interact with Met414; Leu446 and the third residue might be Tyr415 or Val405. From this pharmacophore models, it is expected to have different bioactive conformations rather than Thumb pocket-I, -II. From docking, most compounds are bent on each other to be accommodated in the Palm pocket-I. This leads to more hydrogen bonds and π-π stacking than Thumb pocket-II, especially for compounds having more than one aromatic ring such as compounds 1 and 21-25. This explains the activity of such compounds.

3D-QSAR
CoMFA model calculates steric (S) and electrostatic fields (E) with contributions 0.43, 0.57, while CoMSIA model computes steric, electrostatic, hydrogen bond donor (D), hydrogen bond acceptor (A) and hydrophobic (H) with contributions of 0.10, 0.26, 0.25,0.14 and 0.24, respectively, as given in Table 2. The predictive power of the 3-D CoMFA and CoMSIA models were determined from calculating pIC 50 of the investigated inhibitors. A maximum difference between the actual and predicted of no more than 0.3 indicates excellent models as listed Table 3 [59].
The contour maps, obtained from CoMFA model, Figure 6, show a green region at C7 indicating the requirement for bulky substituents at this region for enhancing potency of the 2 inhibitor. Bulky substitution is essential when green polygon exists and prohibited when it is yellow [57]. The qualitative SAR for compounds 2-6 indicates that the introduction of cyclopentane ring improves the inhibition efficiency compared to cyclobutane, Table 1. However, the sterically disfavored yellow contours for compounds 3, 4, 6 at C16, 17, 19 atoms indicate reduction of the activity upon substitution with bulky groups at these positions. As depicted in Figure 6, for compound 2, the green polygon on the side of the phenol ring will interact with the hydrophobic residue (Leu419, Met423) which increase binding between the inhibitor and receptor. An inspection of Figs. 5-7 indicates that substitution with electron rich groups is expected to improve the inhibitory effect as shown by the red color in Figure 6 for compounds 3,6=>50,52 µM compared to 2,4= >100, 93 µM, Table 1), hydrophobicity increased in Figs. 5,7. This might interpret the increase in the potency when adding aromatic or hetero atoms such as sulfur atom, Table 1. For CoMSIA field, any hydrophobic, acceptor fields are included in the steric field. However, the donor field is embedded in the electrostatic field as concluded from the tripos bookshelf for sybyl (sybyl's manual). The hydrophobic disfavored (gray) over methyl C7 for 2; C17 (gray) disfavored hydrophobic, C4 (orange) favored hydrophobic, C7, C16 (yellow) disfavored steric for inhibitor 3; yellow, gray contours over cyclobutane for compound 4 and finally (yellow) over part of cyclohexane, (orange) favored hydrophobic on C15 for inhibitor 6. These findings may explain comparable activity of compounds 2, 4: of 100, 93 µM and compounds 3, 6: 50, 52 µM as displayed in Figure 4 and Table 1.     Figure 7 shows Steric, electrostatic, hydrophobic, acceptor and donor contour maps around compound 2, 3, 4 and 6, respectively. Sterically favored areas are given in green, sterically unfavored areas in yellow; positive-charge-favored areas in blue, positive-charge-unfavored areas in red. Hydrophobic favored areas given in orange, hydrophobic unfavored areas in gray, donor-favored areas in cyan, donorunfavored areas in purple, acceptor-favored areas in magenta, acceptor unfavored areas in deep red. The maps generated depict regions having scaled coefficients >80% (favored) or <20% (disfavored).

Conclusions
The present work provides an analysis of the interaction between HCV NS5b and 25 compounds for better understanding of the reported potency of these compounds. Molecular docking and 3D-QSAR models were used to study binding of the selected inhibitors with HCVNS5B allosteric sites. Hydrophobic interactions were found to be the most prominent factor for the interaction between the current inhibitors and HCV NS5B. This finding agrees with experiment. Substitutions at particular positions were explored based on electrostatic, steric, and hydrophobic interactions.