Evaluation of interactions between the hepatitis C virus NS3/4A and sulfonamidobenzamide based molecules using molecular docking, molecular dynamics simulations and binding free energy calculations

The Hepatitis C Virus (HCV) NS3/4A is an attractive target for the treatment of Hepatitis C infection. Herein, we present an investigation of HCV NS3/4A inhibitors based on a sulfonamidobenzamide scaffold. Inhibitor interactions with HCV NS3/4A were explored by molecular docking, molecular dynamics simulations, and MM/PBSA binding free energy calculations. All of the inhibitors adopt similar molecular docking poses in the catalytic site of the protease that are stabilized by hydrogen bond interactions with G137 and the catalytic S139, which are known to be important for potency and binding stability. The quantitative assessments of binding free energies from MM/PBSA correlate well with the experimental results, with a high coefficient of determination, R2 of 0.92. Binding free energy decomposition analyses elucidate the different contributions of Q41, F43, H57, R109, K136, G137, S138, S139, A156, M485, and Q526 in binding different inhibitors. The importance of these sidechain contributions was further confirmed by computational alanine scanning mutagenesis. In addition, the sidechains of K136 and S139 show crucial but distinct contributions to inhibitor binding with HCV NS3/4A. The structural basis of the potency has been elucidated, demonstrating the importance of the R155 sidechain conformation. This extensive exploration of binding energies and interactions between these compounds and HCV NS3/4A at the atomic level should benefit future antiviral drug design.


Introduction
Hepatitis C virus (HCV) infection presents a global health threat. Approximately, 200 million individuals are chronically infected with HCV, and the prevalence can be much higher in specific populations [1,2]. It is a major cause of many life-threatening hepatic conditions, including hepatocellular carcinoma, chronic hepatitis C and liver cirrhosis [3][4][5]. Based on the differences in genome sequences, HCV has been divided into seven major genotypes (GT1-7), which are further classified into various subtypes [6].
HCV possesses a single-stranded, positive sense RNA genome that encodes for a single long polyprotein that undergoes proteolytic cleavage to form ten individual proteins: structural proteins (C, E1, and E2), a small ion channel protein, p7, and six nonstructural (NS) proteins (NS2, NS3, NS4A, NS4B, NS5A, and NS5B) [7]. Among these, Jinhong Ren and Tasneem M. Vaid have contributed equally to this work.
there are two proteases, the NS2 cysteine and NS3 serine proteases. The NS2 protease is responsible for cleavage of one site between NS2 and NS3 (shown in a green arrow in Fig. 1a), while the NS3 protease cleaves the four sites downstream between NS3/4A, NS4A/4B, NS4B/NS5A, NS5A/ NS5B [8,9], with the NS4A acting as a critical cofactor in this process [9,10] (Fig. 1a). NS3 is a bifunctional enzyme containing an N-terminal serine protease domain with residues S139, H57, and D81 forming the catalytic triad and a C-terminal helicase domain [11] (Fig. 1b). This NS3/4A protease complex is a key component of the viral replication complex. Along with its involvement in polyprotein proteolysis, it participates in innate immune evasion [12], playing a key role in the perseverance and pathogenesis of hepatitis C in host cells [13].
These NS3/4A protease inhibitors are also used in Direct-acting antiviral agent (DAA) combinations of compounds that individually inhibit the HCV NS3/4A protease, the NS5A polymerase, and/or the NS5B polymerase, and are currently being used to effectively treat hepatitis C [29] in a variety of combinations [30]. Although effective, the therapy has its own shortcomings; (i) The cost of DAAs is quite high, adding to the economic burden on patients [31], and there is at least one report of perceived stress in HCV infected patient under DAA therapy [32], (ii) there have been reports of associated side effects due to off-target engagement [33,34], (iii) not all HCV genotypic variants are equally responsive to the current therapy [35], (iv) resistance-related variants could arise during DAA monotherapy as a result of the high HCV mutation rate [36,37], (iv) failed therapies have been reported for some difficult-to-treat patients, particularly with genotype 3 or cirrhosis complications [38], and (v) recurrences of hepatocellular carcinoma have been reported after HCV therapy with DAAs [4].
Hence, there is a continuing need for discovery of novel and potent small-molecule protease inhibitors of HCV NS3/4A. Most of the reported NS3/4A inhibitors are peptidomimetics with functional C-terminal carboxylate, tetrazole, and acyl sulfonamide substituents [16,[39][40][41]. We have been exploring small molecule functionalities as HCV inhibitor scaffolds [42][43][44] and have found the sulfonamidobenzamide to be a promising scaffold for genotype 1b [43,44]. This study explores computational analysis of the interactions between sulfonamidobenzamide inhibitors and the NS3/4A complex of HCV (genotype 1b) at an atomic level and provides analysis of their binding energies as well as the structural basis of potency differences between them.
Computational modeling of protein-ligand interactions is widely used in modern drug discovery. Molecular dynamics (MD) simulations and binding free energy calculations are efficient tools to investigate thermodynamic properties of proteins and protein-ligand interactions [45] and have been widely used to study various classes of HCV inhibitors [46][47][48][49]. Although binding free energies calculated with the MM/PBSA approach generally do not reproduce the absolute experimental values, they are effective in ranking congeneric inhibitors, and are commonly used to discriminate between weak and strong inhibitors [50][51][52].

Molecular docking of the sulfonamidobenzamide compounds to HCV NS3/4A
A total of twelve compounds with the sulfonamidobenzamide scaffold were selected for comparison, with differing substituents at positions A and C (Fig. 2). All compounds (Fig. 2) were subjected to molecular docking to predict the possible binding poses and to evaluate the binding affinities. The docking protocol was benchmarked by successfully redocking the crystallized ligand F9K [53] into the active site with an RMSD of 0.4 Å (Fig. S1). The docking pose with the highest ChemPLP [54] score for each compound was extracted for comparison. The results show that all inhibitors shared a similar binding pose, with certain differences for aromatic ring A and sulfonamide moieties in the catalytic sites of HCV NS3/4A protease. Ring A tends to occupy two different pockets, S1 or S1', separated by the sidechain of Q41 (Fig. 3a). For example, ring A of F0325-0077, F0325-0086, F0325-0092, F0325-0093, F0325-0125, F0725-0019, F0816-0111, and F1813-0711 occupies the S1 pocket, and can form face-to-face stacking interactions with the amide sidechain of Q41, and edge-to-face stacking interactions with the imidazole ring of H57 and the phenyl ring of F43 (Fig. 3b). However, only edge-to-face stacking with the amide sidechain of Q41 was detected if ring A flipped to the S1' pocket for the binding of F1822-0567, F2708-0139, F2708-0140 and F2730-0247 (Fig. 3b). The sulfonamide moiety forms hydrogen bond interactions with G137 and S139, but different hydrogen bond networks were observed for inhibitors with ring A occupied in the S1 or S1' pocket. For example, for F0325-0077, F0325-0086, F0325-0092, F0325-0093, F0325-0125, F0725-0019, F0816-0111, and F1813-0711, with ring A occupying in S1 pocket, both sulfonamide oxygens are involved in hydrogen bonds with the sidechain OH of S139 and the backbone NH of residue G137. However, this sulfonamide moiety slightly rotates away accompanied by the occupancy of the S1' pocket for ring A in the less active compounds F1822-0567, F2708-0139, F2708-0140 and F2730-0247, and only one oxygen retains the hydrogen bond interactions with G137 and S139. Therefore, the sulfonamide functional group seems important for the binding of these compounds, which is also shown in the binding of the crystallographic ligand F9K (Fig. S2). The connecting phenyl ring B occupies the S3 pocket and remains buried into the hydrophobic pocket for all the compounds' binding predictions, with the binding driven by van der Waal interactions with sidechains of L135, F154, A156, and A157 ( Fig. 3) based on their docking poses. The stacking interactions of the inhibitor amide moiety and aromatic ring C with both sidechains of the catalytic H57 and the helicase domain residue Q526 of HCV NS3/4A were also wellpreserved in all binding predictions. Steric bulky groups were tolerated at ring C, since the broad space of pocket S2 provides space for large groups, allowing for the introduction of more extensive interactions (Fig. 3a). An additional helicase residue, M485, was also involved in van der Waals contacts with ring C in the inhibitors. Therefore, residues M485 and Q526 from the helicase domain of HCV NS3/4A are found to be important for the binding of this inhibitor series from our docking calculations. A hydrogen bond interaction was also found between the backbone carbonyl oxygen of R155 and the amide NH of inhibitors F0325-0077, F0325-0086, F0325-0092, F0325-0093, F0325-0125, F1813-0711,   The chemical structures of twelve sulfonamidobenzamide compounds and their inhibitory activities against the HCV NS3/4A protease. The general scaffold of compounds includes a para-substituted phenyl for ring A (blue), a sulfonamide group, a di-substituted 1,2-phenyl for ring B (green), an amide group, and a substituted thiazole for ring C (red). The IC 50 values for the nine compounds are from our prior work, [a] [44] and [b] [43]. The IC 50 values for the compounds with no superscripts were determined in the current work. Positions 4, 5, 6, 7 are labeled on the R 2 moiety of F0325-0077 (top entry) and F1822-0567. Overall, multiple sulfonamidobenzamides adopt similar binding modes with HCV NS3/4A protease.
Hydrogen bond interactions with G137 and S139 seem more dominant than that with R155, since they were found in all docked inhibitor binding modes. The stacking interactions also played an important role for the binding of these compounds. However, it was still difficult to estimate which interactions predominate in binding, or which interactions were not stable based on docking results alone.

Inhibitory activities against HCV NS3/4A protease
All twelve compounds exhibited various degree of inhibitory activity (IC 50 ) against the HCV NS3/4A protease as summarized in Fig. 2. The IC 50 values of seven compounds were less than 10 µM varying from 1.2 µM to 8.0 µM. Of the seven, F0325-0125 possessed the best inhibitory activity at 1.2 µM. Four compounds (F1822-0567, F2708-0139, F2708-0140, and F2730-0247) showed IC 50 values in the range of 12.7-32.9 µM, while F0325-0077 did not have any activity. Based on these observations, methyl, methoxy and chloro groups are tolerated at the para positions of the phenyl in ring A, while large aromatic substituents are favored for moiety C. Positions 6 and 7 in the benzothiazole ring are more tolerant of large groups than positions 4 and 5, which could be seen from F0328-0092 and F0328-0093. F0725-0019 also shows lower IC 50 values with both positions 6 and 7 in the benzothiazole ring substituted. Methoxy is not S1 S1' S2 S4 S3 Q41 H57 D81 S139 a b Fig. 3 The predicted binding poses from molecular docking. a The surface of HCV NS3/4A is colored in light gray with the substrate binding subsites S4-S1' labeled in red, the catalytic triad in green sticks, and Q41 in orange sticks. The predicted binding poses are shown in overlaid presentation. F0325-0077 (dodger blue), F0325-0086 (grey), F0325-0092 (deep pink), F0325-0093 (cyan), with ring A binding to pocket S1; F1822-0567 (tan), F2708-0139 (sky blue), F2708-0140 (salmon), F2730-0247 (plum) with ring A binding to pocket S1'. b Docking results for twelve inhibitors with HCV NS3/4A. Catalytic triad residues are in green sticks and the other active site residues are in tan sticks. Only the sidechains are shown to depict the interactions clearly, except for residues G137 and R155, for which inhibitors exhibit backbone interactions. Hydrogen bonds are indicated by black dashed lines tolerated at position 4 of the benzothiazole when a methyl was substituted at the para position of the phenyl ring A in F0325-0077. Inhibitory activities vary significantly depending on the sub-structure at the R 2 position.
Furthermore, poor correlation is observed between the ChemPLP score and the experimental values (Fig. 4a), which is a common and major problem with current docking approaches because of protein flexibility and lack of accounting for solvation effects. Thus, MD simulations were applied to offset this disadvantage, to explore the key interactions, and to explain the different activities by calculating the binding free energies using the MM/PBSA approach.

Binding free energy calculations by MM/PBSA approach based on 20-ns MD simulations and correlation with experimental results
Twelve inhibitors were chosen to probe their various interactions with HCV NS3/4A at the atomic level and to explore the correlation between calculated binding free energies and the experimental values. In the MM/PBSA calculations, the explicit water molecules were stripped off, and a continuum implicit water model was used. Interaction entropy was estimated and combined with MM/PBSA computations to obtain more reasonable binding affinities. Our calculations reveal that MM/PBSA showed good correlation (R 2 = 0.92) with experimental data for the 12 inhibitors (Fig. 4b) based on the 20-ns MD simulations and improved obviously comparing with that from docking scores (R 2 = 0.56), which confirm that incorporating protein flexibility as well as solvent effects are essential for protein-inhibitor analysis. In addition, other atomic radii related to the PB [57] model and various GB models, including GB Neck , GB HCT , GB OBC , and GB OBC2 , were systematically investigated in this study, all of which showed good correlation with experimental results (See Supplementary Information Figs. S4 and S5). In general, calculated binding free energies show good correlation with the experimental values.

Hydrogen bond analysis
Hydrogen bonds usually play an important role for the efficacy and specificity in protein-inhibitor interactions. Since docking suggests that several hydrogen bond interactions are involved in these compounds binding with HCV NS3/4A, these interactions were further explored by analyzing the 20-ns MD trajectory of representative compounds to elucidate their stability and contributions. Table 1 lists the three most prominent hydrogen bond interactions for HCV NS3/4A with four representatives, F0325-0125, F0325-0086, F2708-0140, and F0325-0077, possessing different inhibitory activities. Both sulfonamide oxygens were involved in the formation of hydrogen bonds for the most active compounds, F0325-0125 and F0325-0086, with residues G137 and S139 of HCV NS3/4A acting as hydrogen bond donors from docking prediction. These interactions were well preserved with around 90% occupancy throughout the 20-ns MD simulations. Interestingly, the sulfonamide of F2708-0140, showing only one oxygen atom maintaining the hydrogen bond interactions with G137 and S139 from docking, tended to move back with the second oxygen atoms getting involved in hydrogen bond interactions during MD simulations. However, the new joint interaction was not stable, with occupancy much lower than that of F0325-0125 and F0325-0086. For the inactive compound, F0325-0077, despite a network of two oxygen hydrogen bonds with G137 and S139 observed in docking similar to other active compounds, these interactions thoroughly collapsed during MD simulations, which could be the result of an unstable predicted docking pose. The results suggest that the hydrogen bond interactions for the sulfonamide moiety with residues of G137 and S139 are important for maintaining their interactions and inhibitory activity toward HCV NS3/4A. In addition to the hydrogen bond for the sulfonamide core, the hydrogen bond involving the amide linker with the backbone of R155 for F0325-0077, F0325-0086, and F0325-0125 shows only partial occupancy during MD simulations and was predicted not to be pivotal for stabilizing the binding conformations for inhibitors, which can also be seen from the docking prediction that only seven of the inhibitors could form hydrogen bond with R155.

Per-residue energy decomposition to explore the binding hot spots
Per-residue energy decomposition of HCV NS3/4A with representative compounds (F0325-0125, F0325-0086, F2708-0140, and F0325-0077) was investigated to gain further insights into the energetic contributions of individual residues to complex formation using the MM/PBSA (bondi radii) approach. The most favorable interactions are formed by residues in the catalytic sites of the protease and part of the helicase. These are mainly by either hydrophobic interactions or hydrogen bonds. The energy contributions of the more important residues in the HCV NS3/4A-inhibitor interactions are shown in Fig. 5. The decomposition of the enthalpy revealed that all of these complexes have similar interaction networks, which further supports the hypothesis derived from docking that the inhibitors adopt similar binding modes. The most favorable interactions are formed by residues Q41, F43, H57, R109, K136, G137, S138, S139, and A156 in the protease domain as well as M485 and Q526 in the helicase domain. The significant contributions of G137 and S139 indicate their critical role of hydrogen bond interactions with these inhibitors. In addition, residues Q41, F43, H57, and Q526 probe the important function of the stacking interactions. The phenyl ring B occupies the hydrophobic pocket in the HCV NS3/4A protease and was deeply buried in the hydrophobic pocket during the simulations, which can be observed from the A156 contribution. M485 made large contributions to the binding free energies through hydrophobic interactions with ring C. It is notable that residues K136, S139, G137, and F43 contribute much more for F0325-0125 and F0325-0086 than for F2708-0140, and even more than that for F0325-0077. That is consistent with the hydrogen bond analysis stated above for these four compounds during MD simulations with residues G137 and S139, since persistent hydrogen bond interactions are pivotal for their binding. The edge-to-face stacking interaction between the phenyl ring of F43 and ring A is also important to discriminate among the distinct binding activities of the various inhibitors. The interactions and contributions of K136, which exhibits the most favorable binding energy residue for all complexes will be further discussed below.

Alanine scanning mutagenesis for HCV NS3/4A-F0325-0125 complex
Computational alanine scanning was employed to determine the role of active site residues on binding of these compounds. The trajectories of the wild-type HCV NS3/4A binding with the most active compound, F0325-0125, were used to generate the structures of the mutated complex. This Per-residue energy decomposition method depends on the assumption that the effects of mutating a residue to alanine only propagates locally and will result in insignificant entropic changes. For alanine scanning, the above important residues shown in Fig. 5 (except for G137 and A156) were mutated individually, and the results of the mutagenesis are presented in Fig. 6. Changes in the inhibitor-residue interactions associated with the alanine scan show that, in general, mutations of active site residues are highly unfavorable. Significant losses in binding free energies were observed for the catalytic residue of S139, as the S139A mutation removes a key hydrogen bond between the sidechain of S139 and the F0325-0125, responsible for stability and potency. In addition, the mutation of K136A exhibits a substantial destabilization. The hydrophobic and stacking interactions also contribute significantly to the interactions, which can be observed from the energy change by the mutation of Q41, F43, H57, M485 and Q526 to alanine. The mutation of S138A does not alter the binding from this result, since only the backbone is involved in the binding and the sidechain of S138 is away from the binding pocket. Overall, the result suggests that the sidechains of two residues K136 and S139 play pivotal role in the binding of F0325-0125.
Different contributions for K136 and S139 of HCV NS3/4A with F0325-0125 complex K136 and S139 appear to be the two most important residues from per-residue decomposition and their sidechains are indicated to be pivotal for the binding of F0325-0125 to HCV NS3/4A by computational alanine scanning mutagenesis. However, the sidechain roles in nonpolar contributions and polar contributions clearly differ (Fig. 7a). K136 contributes predominantly by van der Waals interactions, which can be distinguished from the snapshot of the most dominant cluster complex from 20-ns MD simulations (Fig. 7b). The K136 sidechain forms tight hydrophobic interactions with ring B of F0325-0125 by its sidechain atoms of C γ and C δ . However, there is no C γ , C δ in the sidechain of K136A, resulting in a dramatic loss of binding affinity for the K136A mutant as well as losing the polar interactions of N ζ with the amide group of F0325-0125. On the other hand, S139 contributes more polar interactions than van der Waals interactions, similar to the other pivotal residue, G137, and dramatically loses ΔΔG for the S139A mutant, which cannot retain this polar hydrogen bond interaction with the sulfonamide oxygen. Therefore, K136A and S139A exhibit the same trend of losing binding affinity, but the losses are derived from different pivotal interactions. Taken together, K136 and S139 show crucial contributions toward vdW and H-bond interactions respectively.

Structural basis for potency differences between F0325-0125 and F0325-0077
The analysis of the 20 ns MD simulation trajectories of F0325-0125 and F0325-0077 docked HCV NS3/4A complexes followed by the comparison of the binding pocket residues of the final frame of 20 ns MD simulation of F0325-0125 with that of the docked receptor complex shows major conformational changes in the side chain orientation of residue R155 (Fig. 8a). For F0325-0125, the flexibility of the thiazole linked 2,5-dichlorophenyl moiety allows the compound to "induce fit" around this flipped orientation (i.e., state II) along with retaining the critical hydrogen bond interactions of the sulfonamide oxygens with the S139 and G137 (Fig. 8b). F0325-0077, which contains the rigid benzothiazole moiety, doesn't fit well around this conformational flip, resulting in a destabilization of the critical H-bonding interactions and as a result, the compound completely collapsed during the 20 ns MD simulation run. Furthermore, the presence of a methoxy at position 4 of the benzothiazole in F0325-0077 results in a potential steric hinderance with the R155 side chain orientation (state II) (Fig. 8c), a conformation found to be stabilized during the 20 ns MD simulation with the F0325-0125 docked complex. The analysis of MD simulation trajectories suggests that the state I conformation of R155 is stabilized by the salt bridge formation between the sidechains of D168 and R155; and hydrogen bonding with the backbone of Q80. (Fig. 8d). The state II conformation is mostly stabilized by the salt bridge formation between the sidechains of D81 and R155 (Fig. 8e).
Comparison with the other crystal structures of HCV NS3/4A complexes also revealed that the R155 residue adopts two different side chain conformations (i.e., orientation around χ 3 ) (Fig. S12). For the macrocyclic ligand bound complexes (PDB IDs 3M5L [58], 6NZT [59], 6NZV [59] and 4A92 [53]), the orientation of the R155 sidechain is different (state I) than that with the PDB structures complexed with non-cyclic compounds (PDB IDs 2FM2 [60], 3LON [61], 3LOX [62]). A closer inspection of the binding pocket reveals that this flipped The key interactions stabilizing the state I and state II conformations of R155 respectively orientation (state II) is most likely triggered by the nearby protonated sidechain of the H57 residue. A protonated histidine side chain group is known to exhibit a high propensity for forming like-charged high-strength contact pairs with protonated arginine [63]. With the crystal structures bound to macrocyclic compounds, the H57 is involved in pi-pi interaction (4A92 [53]) or hydrogen bonding (3M5L [58], 6NZT [59] and 6NZV [59]) with the ligand, which is further stabilized by backbone hydrogen bonding interactions with G137 and R155.

Conclusions
We have analyzed the interactions of a set of sulfonamidobenzamide based inhibitors with the HCV NS3/4A protease using molecular docking, MD simulations, and MM/PBSA calculations. The binding free energies of HCV NS3/4A with the inhibitors using the MM/PBSA method show statistically significant correlation with experimental binding affinities, with R 2 up to 0.92. The hydrogen bond interactions with residues G137 and S139 appears to be pivotal for potency, as compounds F2708-0140 and F0325-0077, which lack stable hydrogen bonding networks at this position, display moderate to no inhibitory activity, with IC 50 values ten to a hundred times lower than those that maintain stable hydrogen bond networks. Residues K136, G137, and S139 make the most favorable contribution to the binding free energy from per-residue energy contribution analysis and were identified as "hot-spots" for binding with the sulfonamidobenzamide based molecules. Additionally, residues Q41, F43, H57, R109, S138, A156, M485, and Q526 made substantial contributions to the binding free energies in all complexes. Furthermore, the sidechains of K136 and S139 respectively show crucial vdW and hydrogen bond interactions with the inhibitors. K136 contributes predominantly van der Waals interactions by its sidechain atoms of C γ and C δ with ring B of F0325-0125, while S139 contributes more polar hydrogen interactions with the sulfonamide oxygen of the F0325-0125. In addition, the structural analysis reveals that the side chain conformation of the R155 residue potentially plays a key role in the efficient binding of these compounds. Thus, this study provides structural and energetic details of the interactions between a set of small sulfonamidobenzamide based inhibitors and HCV NS3/4A, providing a foundation for the design of more potent inhibitors.

HCV NS3/4A inhibitory activity
The full-length genes of wild-type genotype 1b HCV NS3 (HCV polyprotein residues 1027-1657) with a His-tag at the N-terminus and NS4A (residues 1658-1711) with no tag were cloned, co-expressed, and purified as described [64]. All HCV NS3/4A enzyme assays were done in triplicate in black 384-well plates (Matrix Technologies). The HCV NS3/4A enzyme (10 nM final concentration) was prepared in assay buffer (50 mM Tris, pH 7.6, 0.5% Chaps, 15% glycerol, 2 mM GSH, and 0.1 mg/mL BSA), and 20 µL each was dispensed into wells. A series of increasing compound concentrations (0 to 200 µM final concentration at twofold serial dilution) in 100% DMSO were prepared in another 384-well plate. 0.5 µL of varying concentration of compounds were added to 20 µL enzyme solution and incubated for 10 min. The enzyme reaction was initiated by adding 5 µL of a FRET-based substrate Ac-DE-Dap(QXL520)-EE-Abu-ψ-[COO]AS-C(5FAMsp)-NH 2 (Anaspec) at 1 µM final concentration, and enzyme activity was continuously monitored for at least 6 min using a POLARstar OPTIMA microplate reader (BMG LABTECH). The IC 50 values were calculated by fitting the data with the Hill Eq. (1), with SigmaPlot v12, where y is percent inhibition, x is inhibitor concentration, n is the slope of the concentration-response curve (Hill slope), and V max is maximal inhibition from two to four independent assays.

Molecular docking
The X-ray crystal structure of HCV NS3/4A in complex with a macrocyclic protease inhibitor with resolution of 2.73 Å was downloaded from the RCSB protein data bank (PDB entry: 4A92 [53]; unit A was used). For the protein structure, the Protein Preparation Wizard [65] in Schrödinger 2016 was used to prepare the NS3/4A structure for docking, by removing crystallographic waters and ions, filling missing side chains and loops, fixing bond orders, adding hydrogens, assigning partial charges with the OPLS3 force field [66], and minimizing the added hydrogens. The 3D structures of the small molecules were processed using the LigPrep module (Schrödinger 2016) [67]. For each compound, the ionized states and tautomers were generated by using Epik at pH 7.4. All 12 inhibitors (Fig. 2) were docked into the catalytic sites of HCV NS3/4A with F9K removed, and with a radius of 10 Å centered on the catalytic residue of S139 using GOLD v5.2.2 [55]. A rigid receptor and flexible ligands were used during the docking. Standard default settings were used for other parameters. The docked conformation of each ligand was evaluated using the ChemPLP scoring function [54]. The binding pose with the highest score for each molecule was chosen for the following MD simulations. Figures for docking interactions were produced using Chimera [68] and ICM-Molsoft [69].
(1) y = V max x n IC 50 n + x n

MD simulations and MM/PBSA calculations
A total of twelve HCV NS3/4A-inhibitor docking complexes were subjected to MD simulations which were carried out using the PMEMD module of the program AMBER14 [70] to do the energy minimizations and the MD simulations. The RESP atomic partial charges for each of the sulfonamidobenzamide based compounds were generated with the Gaussian09 program at the HF/6-31G* level in the R.E.D. online server [71]. The ff14SB force field was applied to the protein and the general AMBER force field (GAFF) [72] was assigned to inhibitors. Hydrogen atoms for HCV NS3/4A were added with TLEAP in Antechamber. Each system was solvated with a truncated octahedral box of transferable intermolecular potential three-point (TIP3P) explicit water molecules, with at least a 12 Å boundary between the solute atoms and the borders of the box. The appropriate number of Cl − counterions were added to electrostatically neutralize the system. The Particle mesh Ewald method [73] was used for nonbonded interactions with a 10 Å cutoff radius. During the simulation, the SHAKE algorithm [74] was applied to constrain all covalent bonds to hydrogen atoms.
Each system was initially relaxed using energy minimization by 10,000 steps of steepest descent followed with the 10,000 steps of conjugate gradient minimization to remove any steric clashes. This was followed by heating the system from 0 to 300 K with restraints (force constant of 10 kcal mol −1 Å −2 ) on the complex over a period of 100 ps while allowing water molecules and ions to move freely. Subsequently, each system was equilibrated at 300 K and 1 bar for another 100 ps without any restraints. Finally, unrestrained production MD simulations of 20 ns were performed at a temperature of 300 K and a pressure of 1 atm in an NPT ensemble. The coordinates were saved every 5 ps, and a time step of 2 fs was applied. Simulation trajectories were analyzed by cpptraj [73]. MM/PBSA was used to investigate the binding free energies for the complexes from molecular docking. A total of 8,000 snapshots were used to calculate the MM/PBSA. Binding free energies were calculated based on the assumption that the bound and unbound conformations of the protein and inhibitor are quite similar.
The binding free energy between protein and ligand of every snapshot was estimated via Eq. 2.
The free energies of the complex, the protein, and the ligand are characterized as G complex , G protein and G ligand , respectively. The binding free energy (ΔG bind ) for each of the term was calculated using Eq. 3: (2) ΔG bind = G complex − G protein + G ligand (3) ΔG bind = ΔE MM + ΔG sol − TΔS where ΔE MM is the difference in the molecular mechanics energy upon complexation in the gas phase and ΔG sol is the difference in the solvation free energy and −TΔS is the entropy term [75,76].
The molecular mechanics free energy (ΔE MM ) can be further split into van der Waal (ΔE vdW ) and electrostatic (ΔE ele ) energies in Eq. 4: The solvation free energy ΔG sol arises from polar solvation free energy (ΔG pol ) and nonpolar solvation free energy (ΔG nonpol ) as in Eq. 5: ΔG pol represents the polar contribution to solvation and can be estimated by the linearized Poisson-Boltzmann equation using bondi atomic radii and a solvent probe radius of 1.4 Å. The PB calculation used a grid spacing of 0.5 Å, with dielectric constants of 1 and 80 for solute and solvent, respectively. ΔG nonpol is the nonpolar solvation contribution which was determined using an empirical solvent accessible surface area (SASA)-dependent term as in Eq. 6: where γ is the surface tension proportionality constant and β is the offset value.
The entropy term (−TΔS) was computed under interaction entropy calculation method proposed by Duan et al. based on the fluctuations of protein-ligand interaction energies [76][77][78] and is calculated on the 8000 frames extracted from the 20-ns MD simulations by the following Eq. 7:' where β is the reciprocal of KT, and ΔE int pl is the fluctuation of protein-ligand interaction energy around the average energy. Angle brackets represent an ensemble average.

Per-residue decomposition and alanine scanning mutagenesis
Residues that contribute the most to the interactions with the F0325-0125, or the "hot spot" residues, were identified by breaking the effective binding energies into contributions from the individual residues of HCV NS3/4A. The polar contribution to the solvation free energy was determined by applying the MM/PBSA method, using bondi radii.
Then, the specific residue which shows primary contribution from energy decomposition was replaced by an alanine residue by deleting atoms and truncating the mutated residue and calculating the energy change for by this replacement. The MM/PBSA bondi radii were used to calculate the (4) ΔE MM = ΔE vdW + ΔE ele (5) ΔG sol = ΔG pol + ΔG nonpol (6) ΔG nonpol = γ × SASA + β (7) TΔS = KTln⟨e ΔE int pl ⟩ binding free energy differences (ΔΔG) upon alanine mutation. The ΔΔG is the difference between the mutant and wild type complexes defined as Eq. 8, while ΔG mut is the binding free energy for the mutant HCV NS3/4A with F0325-0125 and ΔG wt is the binding free energy for the wildtype HCV NS3/4A with F0325-0125.