Pharmacophore-Assisted Covalent Docking Identifies a Potential Covalent Inhibitor for Drug-Resistant Genotype 3 Variants of Hepatitis C Viral NS3/4A Serine Protease

1. Introduction

Hepatitis C Virus (HCV) is an enveloped, single-stranded, positive-sense RNA virus of the Flaviviridae family [1]. HCV is the causative agent of liver hepatitis C infection, which has an estimated global prevalence of 2.5% [2]. Chronic HCV infections often lead to liver cirrhosis, hepatocellular carcinomas, and liver failure [2,3,4,5]. HCV isolates exhibit vast genetic heterogeneity [6], based on which the viral variants have been grouped into seven genotypes (G1-G7) and various subtypes [2,7,8,9]. G1 and G3 together account for over 80% of the global infections [10,11]. In particular, G1 has a high prevalence in Europe and North and South America [12] and accounts for 46% of the overall HCV infections [11]. Whereas G3 is more prevalent in Australia and South Asia [12] and is causative of 30% of the global HCV burden [10,11]. Other HCV genotypes (G2 - G7) account for the remaining HCV infections and generally remain highly local concerns [11,13].

The HCV genome is 9.6 kb in length, translating into a precursor polyprotein of ~3000 residues [14]. This polyprotein is then enzymatically processed to yield four structural (core, envelope glycoproteins E1 and E2, and an ion channel p7) and six non-structural (NS) proteins (NS2, NS3, NS4A, NS4B, NS5A, and NS5B) [15]. The structural proteins are processed by host peptidases, while the non-structural proteins are cleaved by viral proteases, NS2 and NS3 [16]. In particular, NS2 is an auto-protease that is responsible for intramolecular cleavage between NS2 and NS3 of HCV polyprotein [17,18]. In the case of NS3, the protein has an N-terminal serine protease domain alongside a C-terminal helicase domain with NTPase activity [19,20,21,22]. NS3, together with its cofactor NS4A, makes a heterodimeric NS3/4A complex that cleaves four scissile peptide bonds (3–4A, 4A–4B, 4B–5A, and 5A–5B) between non-structural proteins [23,24,25]. The multifunctional protein NS3 additionally hydrolyses two signaling proteins, including mitochondrial antiviral signaling protein (MAVS) and C-terminal Toll-interleukin-1 receptor (TIR) domain-containing adaptor inducing IFN-b (TRIF) [26]. Importantly, MAVS and TRIF are responsible for activating the immune response to viral infections [26,27]. NS3/4A-mediated processing of viral and host proteins is central to viral replication and innate immune response, rendering it an ideal therapeutic target [19,20,21,22,28,29].

To date, several NS3/4A targeting direct-acting antivirals (DAAs), including boceprevir [30], telaprevir [31,32], narlaprevir [33], simeprevir, grazoprevir, paritaprevir and glecaprevir have been approved by Food and Drug Administration (FDA) [24,25,34] that act to inhibit viral replication [15,35,36,37,38,39,40]. Three of these marketed drugs are covalent (serine-trap) inhibitors [34]. Nucleophilic attack by the hydroxyl group of Ser139 at NS3 protease catalytic triad on the electrophilic warhead leads to the formation of a stable covalent bond between NS3/4A and inhibitor [33,41,42,43,44]. The nature of the electrophile dictates if covalent inhibition is reversible or irreversible [45]. Covalent NS3/4A inhibitors have demonstrated a good safety profile in addition to being potent in their antiviral activity [46,47]. However, the genetic sequence heterogeneity within NS3/4A impacts the efficiency of NS3/4A inhibitors [24,25]. Consequently, the antiviral response varies significantly between HCV genotypes [48]. NS3/4A of HCV G1 exhibits the highest sustained virological response (SVR) rate (up to 99%) to third-generation NS3/4A inhibitors [19,49,50,51,52]. However, the efficacy of DAAs targeting NS3/4A has been shown to be lower in other genotypes [53,54,55,56,57,58]. Specifically, in G3, a predominant HCV genotype in South Asia, including Pakistan [59] – a low to middle-income country, DAA targeting of NS3/4A has presented a major challenge towards affordable HCV treatment [60,61,62,63]. In addition to decreased SVR rates in response to DAAs [64], NS3/4A G3 is associated with an increased risk of liver cirrhosis and hepatocellular cancer [65].

Besides genotypic variation, polymorphism in NS3/4A is also associated with conferring resistance to approved DAAs [66,67]. These single-site mutations [68] within the ligand-binding pocket of NS3/4A G3 [68,69,70,71] cause resistance in drug binding [66,68,71] leading to the emergence of drug-resistant variants [2,9,72]. Nearly 50% of DAA treatment failures have reported resistance-associated substitutions (RASs) [73,74]. The frequently observed NS3/4A RASs include V36, T54, V55, Q80, R155, A156, and D168 [75]. Specifically, A156F/N/V, V36A+R155K/T, V36M+R155T, V36A/M+A156T, T54A+A156S, T54S+A156S/T, and V36M+T54S+R155K have reported resistance to telaprevir [76]. Therefore, the occurrence of RASs within NS3/4A has challenged the clinical success of approved DAAs [77]. Specifically, D168Q polymorphism has rendered NS3/4A G3 “naturally resistant” to most DAAs [67,69,78]. Moreover, the lack of availability of an NS3/4A G3 crystal structure [79,80] further limits the investigation of genotype-specific mutations [68,81] and their impact on NS3/4A’s interaction with inhibitors [82]. Taken together, NS3/4A G3 remains a difficult-to-treat strain.

In recent decades, in silico drug discovery has revolutionised conventional pharmaceutical development [83,84,85]. Various in silico methods assisted in the discovery of potent antiviral compounds against viruses such as dengue [86], zika [87,88], influenza [89], SARS-coronavirus [90,91] and Ebola [92]. Following this, we have combined in silico techniques [91,93,94,95] to identify small molecules targeting HCV NS3/4A G3. In this work, we have modeled NS3/4A variants by considering G3-specific mutations at fourteen positions within ligand-binding pocket, including residues of the catalytic triad (F43L, H57R, Q80K, R123T, I132L, Y134C/R, S139P, R155G, A156T, V158A, C159V, D168Q, C525W/Y, and Q526H/R) followed by pharmacophore-based virtual screening (PBVS) and covalent docking protocol to identify potential covalent inhibitors. Binding interactions of potential antivirals identified in this study were validated using molecular dynamics (MD) simulations. This study reports a potential anti-HCV compound targeting NS3/4A G3-specific mutations exhibiting strong binding affinity scores (-6.5 to -4.1 kcal/mol) that has the potential to be directly tested for in vitro and in vivo studies for HCV drug development.

2. Materials and Methods

2.1. Sequence Retrieval and Modeling HCV NS3/4A G3

In silico modeling was initiated by retrieving sequences of HCV NS3/4A genotype 3 (G3) from the National Center for Biotechnology Information (NCBI) protein database [96] (see Supplementary Data S1). Inhibitor-bound HCV NS3/4A complexes were retrieved from Protein Data Bank (PDB) [97] for binding site analysis (see Supplementary Table S1). The binding pocket of HCV NS3/4A was investigated using LigPlot+ v.2.1 [98] and Maestro v.11 (Schrödinger, LLC, New York) [99], and the ligand-binding residues were identified (see Supplementary Table S2). Multiple sequence alignment (MSA) was performed using Clustal Omega [100], and the identified interacting residues identified were analysed for mutations within G3-specific sequences retrieved from NCBI protein database [96] (see Supplementary Data S2). A representative sequence for each mutation was selected as a G3-specific variant for modeling (see Supplementary Table S3). ESPript 3.0 was used for MSA visualization [101]. Templates for homology modeling were identified from PDB using Protein Basic Local Alignment Search Tool (BLASTP) algorithm [102] (see Supplementary Table S4). Co-crystallized inhibitor-bound HCV NS3/4A complex (PDB ID: 4A92) was selected as a template. Homology modeling of HCV NS3/4A G3-specific variants was performed using SWISS-MODEL Server [103], an automated homology modeling resource.

2.2. 3D Structure Refinement and Validation

The predicted homology models were refined using GalaxyRefine [104]. Quality assessment of the models was done by using FATCAT [105], MATRAS [106], ERRAT [107], ProSA [108], SAVES Server [109], MolProbity [110], Verify 3D [111], PROCHECK [112] and Maestro [99]. Models with the least number of residues in the disallowed region were selected. The root mean square deviation (RMSD) values were calculated using CHIMERA [113] (see Supplementary Table S5).

2.3. Ligand-Based Pharmacophore Modeling and Virtual Screening

Ligand based pharmacophore modeling was carried out with the Phase module of Maestro [114,115] against HCV NS3/4A. A set of over 100 HCV NS3/4A covalent and non-covalent inhibitors, including marketed drugs reported in literature were retrieved from Protein Data Bank (PDB) (see Supplementary Table S6). The features from these ligands were extracted to generate the ligand-based pharmacophore hypotheses. The compounds were screened from small molecule databases including ZINC [116], MolPort [117], PubChem [118], Mcule [119], and ChEMBL [120]. The screening was performed using the Pharmit web server [121] which uses a pharmacophore and molecular shape search. The resulting hits were refined by applying a series of filters, including Lipinski’s rule (RO5), with at least 2 violations [122]. The oral bioavailability assessment was done using Osiris DataWarrior software [123] and pharmacokinetic and pharmacodynamic properties of the screened ligands were evaluated. These strict criteria removed a huge dataset, and the remaining compounds were utilised for more sophisticated covalent docking. The detailed framework applied in the study for targeting HCV NS3/4A G3 variants using in silico approaches is represented in Figure 1.

2.4. Preparation of HCV NS3/4A and Covalent Docking Protocol

The Protein Preparation Wizard in Maestro [124] was used to prepare modeled HCV NS3/4A G3 variants (G3.v1 - G3.v14). The protein structure was pre-processed by assigning bond orders and adding missing hydrogen atoms using Epik [125] at pH 7.0 ± 2.0. Different tautomeric states were generated for each ligand. Moreover, two alternative stereoisomers were generated per ligand while retaining the specified chiralities. Next, the H-bond network was optimised with PROPKA [126] at pH 7.0. The histidine residues within the active site of serine protease were defined including the histidine of the catalytic triad (His57). Moreover, all water molecules were removed, and a restrained minimization was carried out using the newly optimized OPLS3e force field [99],and the convergence criterion of 0.3 Å RMSD was set for all heavy atoms. The compounds resulting from ligand-based pharmacophore search were processed with the LigPrep module of Maestro [127].

Since the active site of HCV NS3/4A contains a catalytic serine (Ser139), it is possible to target it with compounds that bind covalently. For example, within the crystal structure of HCV NS3/4A protease (PDB ID: 3LON), a ketoamide inhibitor narlaprevir is covalently attached to Ser139 of catalytic triad. Therefore, we carried out a covalent docking protocol to screen the resultant compounds using Schrödinger’s Covalent Docking (CovDock) tool [128] implemented in the Maestro molecular modelling package [99]. In addition to Michael addition reaction, we also explored other reaction types, including nucleophilic addition to double bond or triple bond.

The center of the covalent docking site was defined as the centroid of the catalytic Ser139 residue, and a cube grid centered on this point was set with an edge length of ≤ 20 Å. The compounds were initially docked using the fast-virtual screening mode of CovDock in Maestro and Michael addition reaction to investigate the formation of covalent interaction with Ser139. The compounds that successfully interacted covalently with Ser139 were again docked using a thorough pose prediction mode in Maestro [99] and a maximum number of 100 poses was selected. The cut-off at 2.5 kcal/mol was set to retain poses for further refinement. The poses were ranked using Prime [129] MMGBSA score that employs variable dielectric generalized Born solvation model (VSGB 2.1) [130,131] in the OPLS3e force field [132]. These steps were repeated for other reaction types including nucleophilic addition to a double and a triple bond, nucleophilic substitution, aryl and nitrile activated conjugate addition to the alkyne. The top-ranking compounds were selected for validation using Molecular Dynamics (MD) Simulations.

For G3.v3 with S139P mutation, molecular (noncovalent) docking-based virtual screening was performed to identify potential protease inhibitors. The variant was prepared using a protein preparation wizard in Maestro [99]. Hydrogens were added, and proper bond order and energy minimization were performed. The molecular interaction grid was centered at the active site containing the catalytic triad and was prepared using the receptor grid generation module in Maestro [99]. Molecular docking was performed on selected compounds from pharmacophore-based virtual screening using Glide program of Schrodinger [133]. The resulting lead compounds were analysed for the binding interaction with the catalytic triad.

2.5. Molecular Dynamics Simulations

Molecular dynamics (MD) simulations of each complex were performed to investigate the binding potential of the hit compounds interacting with the residues inside the active site of HCV NS3/4A G3. MD simulations were performed for 100 ns using AMBER 20 simulation package [134]. Simulation trajectories elucidated complex stability, and interaction profiles were investigated. MD simulation protocol as described previously [135,136,137] was implemented however, the length of the production run was increased from 20 to 100 ns. Parameterisation was performed, and a new molecular topology file for each ligand was created. This was initiated from each ligand and ended in covalently bound Ser139. The topology and coordinate files of the complexes were generated using tleap program of AMBER. The Antechamber package of AmberTools was utilised, and parameters were extracted from the GAFF force field (GAFF) [138]. Counter ions were added around the ligand-protease complex to neutralise the charges of each simulation system. The complex was centred in a dodecahedral TIP3P [139] water box with a distance of 10 Å between the solute and the box edge. To maintain constant bond length, covalent bonds were constrained using the SHAKE algorithm [140]. The system was heated and equilibrated after a stepwise minimisation. A production run was performed at 300 K and 1 bar pressure for a period of 100 ns. The time step of 2 fs was set and the trajectory snapshots were saved every 2 ps for onward analysis using the CPPTRAJ program [141] of AMBER.

2.6. Binding Free Energy Calculations

The binding free energies (ΔGbind) of HCV NS3/4A G3 variants complexed with the most promising hit compounds were calculated using the MM-GBSA method, implemented in AMBER 20. For each system, 10000 snapshots were generated from the last 50 ns stable trajectories with an interval of 5 ps. The total binding free energy was calculated as a sum of solvation free energy (ΔGsol) and the molecular mechanics binding energy (ΔEMM) as given below:

ΔE_MM = ΔE_int + ΔE_ele + ΔE_vdw

ΔG_sol = ΔG_pl + ΔG_np

ΔG_total = ΔE_MM + ΔG_sol

ΔG_bind = ΔE_MM + ΔG_sol – TΔS

Here, ΔEMM includes electrostatic energy (ΔEele), internal energy (ΔEint), van der Waals energy (ΔEvdw), and the polar (ΔGp) and non-polar (ΔGnp) energy components contributing towards total solvation free energy (ΔGsol). ΔGtotal is the free energy of binding evaluated for both MM-GBSA and MM-PBSA methods after entropic calculations (-TΔS). Per-residue energy decomposition analysis was performed using the MM-GBSA method to estimate the contribution of interacting residues towards ligand-binding. The binding energy was calculated as ΔGresidue using the equation below:

ΔG_residue = ΔE_MM + ΔG_sol

Here, the ΔGresidue denotes the total energy obtained from sidechain and backbone energy decomposition. Only those amino acids were considered within 8 Å of the active site. The sum of energy contributions from each residue is equal to the system’s overall binding energy [142]. This strategy has been well established for protein–ligand binding studies and for studying protein–protein interactions [95,143].

3. Results and Discussion

3.1. Analysis of Sequence Variation and Molecular Modeling

The results from sequence alignment revealed the catalytic triad of HCV NS3/4A G3 as a site of viral polymorphism due to mutations including H57R and S139P (Figure 2). The analysis further showed mutations within HCV NS3/4A ligand-interacting residues at positions 43, 80, 123, 132, 134, 139, 155, 156, 158, 159, 168, 525 and 526. Of these residues, mutations at positions 80 [144,145,146], 155 [68,147], 156 [148], 123 [66,67] and 168 [149,150,151,152] have been frequently reported as drug resistance-associated substitutions (RASs) in HCV NS3/4A. A representative sequence for each mutation was selected to model the HCV NS3/4A G3 variant (see Supplementary Table S3).

Crystal structure of HCV NS3/4A protease-helicase G1b (PDB ID: 4A92) with the sequence identity in the range 74% - 82% (E-value 0.0) for 15 HCV NS3/4A G3 variants was selected as a template for homology modeling (see Supplementary Table S4). The structure was minimized and used for the HCV3/4A modelling (see Supplementary Figure S1).

The predicted structure quality score was computed using ERRAT [107] and Verify 3D [111] and ranged between 91.33 to 94.85. The MolProbity [110] quality values ranged between 1.50 - 1.61. The superimposition of variants with the template indicated root mean square deviation (RMSD) values in the range 0.28 Å - 0.30 Å (see Supplementary Table S5). The Ramachandran plots generated for each variant model showed model residues >95% in the favored regions (see Supplementary Figure S2). Furthermore, we evaluated the reliability of HCV NS3/4A G3 wildtype (WT) model with its respective templates (PDB IDs; 4B75, 1CU1, 3O8B, 5WDX and 4B6E) and root-mean-square-fluctuations (RMSF) were compared after 100 ns MD simulations. The results indicated no significant fluctuations and the movements of secondary structural elements were analogous to the temaplates, while the catalytic triad remained converged below 1Å (see Supplementary Figure S1).

3.2. Ligand-Based Pharmacophore Search

Next, we incorporated both covalent and non-covalent features of reported HCV inhibitors to build a pharmacophore hypothesis against HCV NS3/4A G3. To identify compounds with the most likely interacting substructures and the potential warheads that can establish covalent bonds with catalytic Ser139, a set of well-defined HCV NS3/4A inhibitors was retrieved (see Supplementary Table S6). Using the PHASE 4.0 module of Schrodinger molecular modeling [99], pharmacophore hypotheses were generated. Five hypotheses were generated, and the hypothesis with the best survival score (5.080) and fitness score was selected. This model consisted of four features (AAHH), including two hydrogen bond acceptors (A4 and A6) and two hydrophobic groups (H14 and H15). The pharmacophore hypothesis was validated with receiver operating characteristic (ROC) analysis [153] to assess its ability to correctly classify compounds as active or inactive [154]. The performance of the pharmacophore hypothesis was further evaluated by the area under the curve (AUC) of the corresponding ROC curve (see Supplementary Figure S3). The validated pharmacophore model AAHH curve (see Supplementary Figure S4) was used to screen databases in the Pharmit web server [121] containing 341,098,760 compounds with 1,606,918,818 conformations. Based upon Lipinski’s rule of five (Ro5), 1,321,302 drug-like compounds were obtained which were further filtered based on ADMET (absorption, distribution, metabolism, excretion, and toxicity) properties in Osiris DataWarrior software [123]. Unique compounds from the filtered 4215 compounds were prepared for docking studies (see Supplementary Table S7).

3.3. Covalent Docking-Based Virtual Screening

Ser139 is a key residue in the active site of HCV NS3/4A, which makes it an attractive target for the design of covalent inhibitors [33,155,156]. Covalent docking-based virtual screening was performed to target the reactive nucleophilic Ser139 to identify compounds with reactive electrophilic moieties. The virtual screening was conducted in two steps, involving i) the selection of candidate ligands based on their relevant conformation such that reactive groups are in close proximity to Ser139; ii) and a virtual chemical reaction between the reactive groups, leading to the formation of a stable covalent bond (S-C) [157,158]. Ligand poses within 5 Å distance cut-off from Ser139 were kept, forming a covalent bond (S-C) according to the reaction type. To identify potential covalent inhibitors from a pharmacophore-based screened library containing chemical warheads, the covalent binding reactions used included Michael addition, nucleophilic addition to a double and a triple bond, nucleophilic substitution, and aryl and nitrile activated conjugate addition to alkyne. The free energies of binding were calculated with the Prime/MM-GBSA method for all the docked poses (see Methods section for details). The compounds with the lowest Prime/MM-GBSA and/or CovDock scores were considered for further molecular inspection.

Overall, a total of 23 compounds from the results of the Michael addition screening, while no compounds could be identified from other reaction types: nucleophilic addition to double bond and nitrile activated conjugate addition to alkyne. Inspection of the structures and binding poses led to the selection of 7 compounds (see Supplementary Table S8 and Supplementary Figure S5). Amongst these, only one inhibitor candidate, CHEMBL569970 (cpd-217: PubChem45485999) (Figure 3) reported high CovDock docking score with all variants incorporating G3-specific mutations and was selected for a detailed interaction analysis. A detailed interaction analysis of cpd-217, and its stability inside the active site of HCV NS3/4A G3 variants, was assessed through 50 ns MD simulations and binding free energies were calculated.

Molecular docking of G3.v1 resulted in 60 compounds (see Supplementary Table S9 and Supplementary Figure S6). Ligand interaction analysis was performed to filter leads. The compound cpd-217 indicated stable interactions, which were validated by molecular dynamics simulation analysis.

3.4. Molecular Insights of Identified Potential Covalent Inhibitor (cpd-217)

The docked complexes of a potential hit with HCV NS3/4A variants were investigated for molecular interactions governed by covalent bond formation. The structures of the lead compound and their molecular interactions with the binding pocket of HCV NS3/4A G3 variants, v1 to v14 (Figure 4A-N), along with the corresponding reaction sites between the chemical warhead of the lead compound and the catalytic Ser139 are displayed in Figure 4.

Overall, cpd-217/ G3 complexes, the catalytic serine (Ser139) was covalently bound to the Michael acceptor warhead (≤ 1.9 Å distance), except HCV NS3/4A G3.v3, which had proline instead of serine at position 139. The residues within the binding pocket [159] of all variants interacting with cpd-217 were compared with the residues interacting with the alpha-ketoamide inhibitor, boceprevir (PDB: 2O8) (see Supplementary Figure S7). Like Boceprevir, cpd-217 formed interactions with Gln41, His57, Ile132, Leu135, Lys136, Gly137, Ser138, Ser139, Phe154, Arg155, Ala156, Ala157, and Val158. Additionally, cpd-217 interacted with Thr42, Phe43, Val55, Gly58, Asp81, Val158, Cys159, Met485, Phe486, Asp487, Ser488, Val524, Cys525, Gln526, and His528 (Table 1). All variants exhibited R123T (except R123S in G3.v10), I132L, Y134C (except G3.v2) and D168Q mutations of which R123T and D168Q are frequently reported as RASs. Additionally, in G3.v1 with H57R mutation, cpd-217 made five H-bonds with Gln41, Thr42, mutated catalytic His at position 57 to Arg, Leu135 and Arg155. In G3.v2 with Q80K mutation, cpd-217 made four H-bonds with Gln41, His57, Leu135 and Gly137. G3.v3 exhibited catalytic Ser139 mutation to Pro, besides R155G mutation. Cpd-217 made three H-bonds with Asp81, Asp487 and Gln526 and established a salt-bridge interaction with Asp81. In G3.v4 (C525W) and G3.v5 (C525Y), cpd-217 made four H-bonds with Gly137, Ser139 and Ala157. Additionally, a π-π interaction was observed with His57 in G3.v4. In G3.v6 with the F43L mutation, cpd-217 made three H-bonds with residues Thr42, Gly137 and Ala157. Cpd-217 was observed making two H-bonds with Leu137 and Gly137 in G3.v7 (Q526H) and an additional H-bond with His57 in G3.v8 (Y134R). In G3.v9 (V158A), cpd-217 made π-cation interaction with Lys136 besides interacting with Gln41, His57 and Arg155 making three H-bonds. Cpd-217 made four H-bonds with His57, Leu135 and Ser139 besides covalently inhibiting Ser139 in G3.v10 with R123S mutation. In G3.v11 with A156T mutation, no H-bond formation was observed. However, cpd-217 made six H-bonds with Gln41, His57, Leu135 and Ser139 in G3.v12 exhibiting Y134S mutation. A salt-bridge interaction between cpd-217 and Asp81 was observed in G3.v13 with Q526R mutation. Additionally, two H-bonds between cpd-217 and Asp81 and Gly137 were observed. Cpd-217 made five H-bonds with Gln41, Gly137, Ser139 and Ala157 in G3.v14 with Y134T and C159V mutations.

These H-bonds were evinced from the reported crystal complexes, where eight hydrogen bonds are conserved between protease residues (Gly137, Ser138, Ser139, Arg155, Ala157, Ser159) and the viral substrate residues [160] indicating that they played a vital role in binding within wild type (WT) systems [34].

We observed that mutations within the catalytic triad residues: His57Arg and Ser139Pro hindered interaction of cpd-217 with Gly137, a residue of the oxyanion hole, in G3.v1 and v3. Mechanisms of drug resistance due to mutations at positions 80 [144,145,146], 155 [68,147], 156 [148], 123 [66,67] and 168 [149,150,151,152] have been well-investigated [161,162] and the importance of these residues in determining the specificity of inhibitors is substantial. RASs R123T and D168Q in addition to mutation at catalytic triad residue H57R (G3.v1), RASs R123T and Q80K (G3.v2) and RASs R123T and R155G along with mutation at catalytic S139P in G3.v3 inhibited interaction of cpd-217 with Gly137. Moreover, R123T and D168Q along with another RAS A156T weakened the association of cpd-217 with G3.v11 resulting in no H-bond formation or hydrophobic interactions (Figure 5A-N).

However, the success of virtual screening is delineated by finding novel chemical structures, new scaffolds are clearly preferred compared to already known scaffolds. Moreover, due to the ability of covalent inhibitor to bind irreversibly to off target proteins, which could lead to toxicological effects, such as immune responses. We utilized SEA (Similarity Ensemble Approach; http://www.sea.bkslab.org) server to identify the presence of similar chemical scaffolds of cpd-217, based on the set-wise chemical similarities against ChEMBL database. The results indicated no structural similarity with any known anti-viral inhibitors (Supplementary Table S10).. We further predicted the unfavorable side-effects due to off-target effects of known molecules and drugs, SwissTargetPrediction web server (http://www.swisstargetprediction.ch) was utilized which combines different measures of chemical similarity based on both chemical structure (2D) and molecular shape (3D). The cpd-217 was found to have less than 0.1 off-target probability based on cross-validation analysis in ChEMBL databse for human protein ligands (Supplementary Table S11).

3.5. The Stability and Flexibility of cpd-217 through MD Simulation

Molecular dynamics studies have emerged as a reliable method for investigating the stability of protein-ligand complexes [163,164,165,166,167]. Recently, these methods have been combined with more sophisticated binding free energy calculations (e.g. MM/GBSA) to explore the drug resistance mechanisms of HCV resulting from several key mutations of NS3/4A [57,58,168,169]. Here, docking energies expounded only the initial conformation of the cpd-217 at the active site of NS3/4A G3 variants. These complexes were further exploited through MD simulations to predict the most likely binding mode between cpd-217 and the associated residues in all G3 variants (v1 to v14). The present study is primarily involved in the conformation of cpd-217 in G3 variants. Therefore, we were interested in analysing the stability and conformational flexibility of the cpd-217 bound to each G3 variant and corresponding energy contributions of potential residues towards ligand. To elucidate the dynamic stability and to ascertain the rationality of the ligand sampling, root-mean-square deviation (RMSD) and root-mean-sqaure-fluctuation (RMSF) values of protein backbone atoms of each G3 variant and heavy atoms of cpd-217 relative to the respective initial structures were calculated, and RMSD/RMSF trajectories were analysed throughout 100 ns. The representative RMSD and RMSF trajectory plots of each G3 variants with bound cpd-217 are displayed in Figure 6A-N and Figure 7A-D.

Overall, the RMSD values of protein backbone atoms (in all variants) and heavy atoms of cpd-217 remained within 3Å throughout the simulation period, except minor variations (slightly higher than 3Å) were observed in G3.v7 and G3.v9 within the initial 30 ns. The RMSD of backbone atoms of the order of 1 – 3 Å with no significant conformational change confirmed that the system is well equilibrated and cpd-217 remained stable inside the binding pocket of all NS3/4A variants during the simulation period [170,171,172]. Whereas, RMSFs highlighted the flexible regions of G3.WT and G3 variants. There were no pronounced Cα-RMSF differences occurred except few small fluctuations for residues 65 to 72 and 85 to 90. Moreover, no gradual fluctuations were observed in the the catalytic site, which evinced the favourable conformation of cpd-217 inside the binding pocket and depicted convergance of catalytic triad residues (< 1 Å) throughout the simulation period.

3.6. Binding Free Energy Calculated by MM/GBSA Method

The MM/GBSA method has been widely employed in improving the protein-ligand docking results [171,172,173] considering the protein backbone dynamics, electrostatic, van der Waal (vdW), and entropic contributions on the overall binding energy of the complex. The MM/GBSA approach has been widely used in expounding the mechanisms of mutation-induced drug resistance [58,174,175,176,177,178,179,180]. MMGBSA total binding free energy values were extracted from initial frames (avg. of 1000 snapshots from 0.01 to 0.5 ns and referred to as “before MD” or “initial”), and final frames (avg. of 10000 snapshots from last 50 ns, after every 5.0 ps and referred as “after MD”). The total binding free energies (ΔGtotal) and energetic components are tabulated in Table 2. We also calculated the affinity score variation (ΔAS) between wild-type (WT) and G3 variants scores to examine the impact of mutation in overall binding with the cpd-217. To further evaluate the binding stability during simulation, energy shift (ΔΔGshift) was incorporated, which supported the stabilising effect of mutation on ligand binding in protein/ligand complex formation, i.e., more negative value represents more favourable stabilising effect in complex with cpd-217. According to the MMGBSA values during MD simulations (Table 2), the overall binding free energies (ΔGbind) of WT was -14.4 kcal/mol. The major energy contributions to the ligand-binding were due to the favorable van der Waals interactions for WT/cpd-217 complex (ΔEvdw = -41.7 kcal/mol) as compared to electrostatic interactions (ΔEele = -32.8 kcal/mol). In contrast, total solvation energies (ΔGsol) showed unfavorable contributions. These values of WT were compared with the binding energies of all G3 variants. The lower binding energy agrees to a higher affinity towards complex stability and vice versa [173,181]. All G3 variants (v1 to v14) possess few common mutations at Arg123, Ile132 and Asp168, along with other specific mutations, the quantitative information of important residue’s contribution towards cpd-217 binding was exploited through per-residue decomposition analysis (Figure 8).

In Table 2, the ΔGbind of the WT complex was more favorable than the G3 variants, except G3.v2 and v14 (-15.6, and -15.1 kcal/mol). On the contrary, the G3 variants with lower binding affinities indicated that the specific mutations in these variants will trigger drug resistance to a certain extent; this phenomenon has been explicitly elaborated in several molecular modeling assisted drug-resistance studies [53,180,182,183,184,185,186]. Likewise in WT, it was obvious that the ΔEvdw contribution is the main component in the total binding affinity in each variant (ranged from -33.5 to -44.2 kcal/mol), and the following was the ΔEele (ranged from -26.4 to -32.5 kcal/mol). The total solvation energies (ΔGsol) that counteract the electrostatic interactions of all systems were within a difference of < 3 kcal/mol from the WT (36.4 kcal/mol). The contributions of the conformational entropy (-TΔS) in all variants ranged between 15.2 and 25.9 kcal/mol. It was observed that the conformational entropic contributions for the complexes have no impact on the order of the ligand’s free energy of binding (ΔGbind). Therefore, the van der Waals contribution was considered more crucial for cpd-217 interaction with NS3/4A protease and differentiating the binding affinities among these variants.

Amongst the G3 mutants, only G3.v2 and G3.v14 displayed a negative affinity variation (ΔAS) of -1.2 and -0.7 kcal/mol, respectively and depicted slightly higher binding affinity with cpd-217. All other G3 variants displayed a positive affinity variation (ΔAS), as low as +2.8 and as high as +8.1 kcal/mol, and exhibited lower binding affinity than WT. Particularly, the ΔGbind of G3.v4 to G3.v7, and G3.v9 to G3.v11 were approximately the same (ranged from -9.1 to 9.9 kcal/mol, with < -1 kcal/mol variation) and exhibited similar positive correlation (ΔAS ranged from +4.6 to +5.3 kcal/mol). The trend was likely due to the presence of four combined mutations (R123T + I132L + Y134C + D168Q) in these variants (except G3.v10 with R123S mutation), which triggered a similar impact on cpd-217 binding as can be observed from converged protein backbone deviation (<3 Å) (Figure 6). Overall, the variants with R123T + I132L + Y134C + D168Q mutations displayed low binding affinity compared to the WT. More in-depth per-residue decomposition analysis of these variants further supported the underlying impact of these mutations upon ligand binding (Figure 8). Among these G3 variants (G3.v4 to G3.v7, G3.v9 to G3.v11), the interaction energy (ΔGresidue) between cpd-217 and mutated residues R123T (-0.99 to -1.84 kcal/mol, while -0.44 kcal/mol in R123S in G3.v10), I132L (-0.95 to -1.46 kcal/mol), Y134C (-0.23 to -0.63 kcal/mol), D168Q (-0.14 to -0.61 kcal/mol) were less favorable than those in the WT complex (-2.14, -1.56, -0.96, -0.95 kcal/mol). Moreover, R123T + I132L + Y134C + D168Q mutations in these variants (v4 to v7, v9 to v11) further reduced the energy contributions by Gln80 (-0.8 to -1.45 kcal), Asp81 (-0.87 to -1.51 kcal/mol), and Ser139 (-1.72 to -2.13 kcal/mol) as compared to its WT complex (-1.89, -1.56, and -2.36 kcal/mol).

G3.v3 displayed the highest positive affinity variation (ΔAS = +8.1 kcal/mol) and showed the lowest binding affinity of -6.3 kcal/mol among all G3 variants. The underlying impact of G3.v3 on cpd-217 was mainly due to the S139P mutation and combined effects of R155G and D168Q mutations. Note that the essential function of Asp168 is to stabilise the conformation of Arg155 to maintain the favorable interaction between Arg155 and ligand [53,54,58,82,162,168,187,188,189]. On the other hand, Proline residues are known to cause short- and long-range disruptive changes in secondary structural elements by causing steric hindrance [190]. Ser139 is one of the important catalytic residues and substitution of this residue has been reported in altered protease activity [191]. An incorporation of glycine at position 155 completely diminished the important interactions with ligand. At the same time, the substitution of S139P further augmented the impact on cpd-217 binding, as evident from the lowest energy shift (ΔΔSshift) of only -2.3 kcal/mol. The contribution energies by these substituted residues further delineated the binding profile with cpd-217. In G3.v3, the interaction energy was abruptly reduced between cpd-217 and mutated residues Pro130 and Gly155 (-0.46 and +0.86 kcal/mol) compared to Ser139 and Arg155 (-2.36 and -2.45 kcal/mol) in WT complex. Moreover, it also impacted the other catalytic triad residues where the interaction energy (ΔGresidue) of His57 and Asp81 were reduced to -1.56 and -1.14 kcal/mol compared to WT complex (-2.8 and -1.56 kcal/mol). Hence, the binding affinity of G3.v3 was the lowest compared to WT. Likewise, H57R in G3.v1 exhibited a similar impact on the catalytic triad, where the incorporation of the long side chain of Arg57 obstructed the binding site which posed the cpd-217 in different conformation as can be seen from a sudden fluctuation at ~50 ns in Figure 6A, that remained stable afterwards. The Arg57 displayed less interaction energy (-1.42 kcal/mol) as compared to His57 in WT, and the conformation of the side chain of Arg57 further reduced the interaction energy of Asp81 (-0.24 kcal/mol) and Ser139 (-0.84 kcal/mol) as compared to WT complex.

G3.v2 and v14 showed negative affinity variation among all G3 variants and displayed slightly higher binding affinities of -15.6 and 15.1 kcal/mol compared to WT. In G3.v2, the Q80K slightly increased the interaction energy of adjacent Asp81, from -1.56 (WT) to -1.8 kcal/mol, while in G3.v14, the Y134T instead of Y134C (as seen in other G3 variants) along with other mutations triggered a combined effect towards cpd-217 (Table 2). The underlying impact in these two variants (G3.v2 and v14) with bound cpd-217 was evident from the highest energy shift (-8.7 and -9.9 kcal/mol, respectively) during the simulation period, which had a stabilising effect.

To summarize, the structural and pharmacophoric characteristics described here can be utilised to find new leads in compound databases and to design new inhibitors that target G3 variants. This work has supported the idea that further optimization of structural and pharmacophore properties can lead in developing a multifunctional small-molecule inhibitor that targets all of the common G3 specific mutations. Future efforts would be necessary to create novel multi-functional anti-HCV inhibitors to study the structural and pharmacophore properties of cpd-217 that are responsible for its diverse activity against all G3 variants.

4. Conclusions

This work reports the drug-resistant mutations within ligand-binding residues including catalytic triad of HCV NS3/4A. The reported mutations are specific for HCV NS3/4A genotype 3 (G3) prevalent in developing countries, including Pakistan. We report mutations within the catalytic triad residues H57R and S139P and RASs Q80K, R123T, R155G, A156T and D168Q in HCV NS3/4A G3 hinder the interaction of drugs with Gly137 of oxyanion hole. This work highlights the pharmacoinformatics approaches utilised to identify a potential covalent inhibitor of HCV NS3/4A G3 to treat Hepatitis C Virus. Several ligands and FDA-approved marketed drugs were used to generate a pharmacophore with a similar scaffold to screen multiple small molecule libraries. Pharmacophore-based virtual screening (PBVS) followed by covalent docking protocol identified cpd-217 (CHEMBL569970; PubChem45485999) as a potential inhibitor of HCV NS3/4A G3 serine protease. The binding affinity, molecular interactions and stability of binding of the lead compound were investigated by molecular docking protocol and molecular dynamics simulations analysis. The potential warhead identified in this work can serve as a guideline to design covalent inhibitors targeting the catalytic Ser139 considering G3-specific drug-resistant mutations within HCV NS3/4A. The proposed inhibitor may play a key role in expediting the drug discovery process and can be tested in clinical trials to treat HCV. This approach can provide a plethora of energetic information, including the binding free energy between protein and ligand, in addition to enriched structural–dynamical information of protein complex structures in solution. Such information is critical for understanding the nature of protein-ligand interactions and guiding drug design and development, which experimental techniques struggle to readily provide.5. Patents