A Comprehensive Computational Analysis for the Binding Modes of

Article pubs.acs.org/journal/aidcbc

A Comprehensive Computational Analysis for the Binding Modes of Hepatitis C Virus NS5A Inhibitors: The Question of Symmetry Marawan Ahmed,† Abhishek Pal,§ Michael Houghton,‡,# and Khaled Barakat*,†,‡,# †

Faculty of Pharmacy and Pharmaceutical Sciences, University of Alberta, 116 Street & 85 Avenue, Edmonton, Alberta, Canada T6G 2R3 § Indian Institute of Technology, Kharagpur, West Bengal 721302, India ‡ Li Ka Shing Institute of Virology, University of Alberta, 116 Street & 85 Avenue, Edmonton, Alberta, Canada T6G 2R3 # Li Ka Shing Applied Virology Institute, University of Alberta, 116 Street & 85 Avenue, Edmonton, Alberta, Canada T6G 2R3 ABSTRACT: Direct-acting antivirals (DAAs) form the current standard of care (SOC) against hepatitis C virus (HCV). These drugs selectively target the viral proteins, offering a unique mechanism to avoid toxicity, to increase their efficacy, and to evolve from decades of interferon- and ribavirin-based therapy. Among the promising HCV targets for DAAs is the NS5A protein, and daclatasvir (DCV) forms a first-in-class compound that selectively targets this protein. Despite the exceptional potency of DCV (∼picomolar IC50) and although several DCV derivatives have been approved for human use or are close to approval, the exact mode of action of these drugs is still incomplete. This is simply due to the vast complexity of cocrystallizing DCV with NS5A in the absence of two amphipathic helices that are required for DCV binding. In this context, computational modeling provides a unique alternative to solve this problem. Here, we build upon our recent discovery of a completely symmetrical interaction between DCV and NS5A and investigate the mode of binding of six other structures similar to DCV. The selected compounds include both symmetric and asymmetric molecules. In addition, we show that our model correlates very well with mutations that can confer resistance to DCV. The current study enhances our understanding of the mode of action of this class of HCV inhibitors and helps in defining the origin of resistance to these drugs. KEYWORDS: HCV, NS5A, direct-acting antivirals, Daclatasvir (DCV), symmetric ligands, asymmetric ligands, NS5A mutations

H

Despite this exceptional potency, we know very little about the exact mode of action of these compounds, making them one of the most enigmatic classes of inhibitors.10 Several mutational analyses mapped the interaction of DCV and its derivatives to the first 100 amino acids of NS5A. These mutations can dramatically reduce the activity of the drug by 10−10000-fold.11 Searching for clues to understand the precise binding of these drugs has been a critical research area since their discovery, and it was only very recently reported that DCV and DCV-like structures can indeed directly bind to the NS5A protein in solution.12,13 Although this discovery ended a huge scientific debate regarding the mode of action of these drugs,14,13 how DCV really interacts with NS5A remains a crucial question that still needs to be answered. It is practically impossible to obtain a cocrystallized NS5A inhibitor with the NS5A protein.12,15 This simply stems from the fact that the two N-terminal amphipathic helices (residues 1−25) are absent in every available crystal structure of NS5A. Given the importance of these helices in DCV-NS5A binding, as described below, a

epatitis C virus (HCV) is still a major global health problem. Chronic infection with HCV can lead to several liver complications, including hepatic cirrhosis and hepatocellular carcinoma. With no HCV vaccine currently available, the standard of care (SOC) treatment of HCV for the past three decades has been limited to interferon and ribavirin therapy.1,2 However, due to their severe side effects and restrictions in treating the different HCV genotypes, the world has been seeking a new type of treatment, a treatment that is free of interferon and ribavirin and targets the viral products, not the host proteins.2,3 A number of these drugs have been already approved, revolutionizing the landscape of HCV treatment and switching the HCV SOC to their advantage.4,5 They have been termed direct-acting antivirals (DAA), and they target mainly three viral proteins, namely, the NS5B RNA polymerase, the NS5A protein, and the NS3 protease.4,5 Moreover, combining these drugs has led to unprecedented synergistic effects and offered, for the first time, a pan-genotype activity, providing an effective cure for all HCV genotypes.6,7 Among these promising drugs is Daclatasvir (DCV; also known as BMS-790052 and Daklinza) and its derivatives. DCV represents a first-in-class NS5A inhibitor with a staggering picomolar inhibition activity and a clinical efficacy of a ∼3−4 log10 HCV RNA reduction in infected patients using only DCV as a monotherapy.8,9 © 2016 American Chemical Society

Special Issue: Host-Pathogen Interactions Received: June 20, 2016 Published: August 3, 2016 872

DOI: 10.1021/acsinfecdis.6b00113 ACS Infect. Dis. 2016, 2, 872−881

ACS Infectious Diseases

Article

HCV inhibitors and helps in defining the mechanism of resistance to these drugs.

realistic structural study must incorporate these residues to truly understand this interaction. With this big experimental challenge, computational modeling offered a new insight into overcoming this hurdle. Here, we continue our efforts toward understanding how DCV and DCV-like structures interact with NS5A. In our previous modeling study, we described a model for NS5A− DCV binding with the main conclusion that DCV binds only to the closed “back-to-back” NS5A dimer and that it has to bind in a symmetrical mode of binding.16 Whereas a number of studies supported DCV binding to the closed NS5A conformation, other modeling studies did not exclude the possibility of DCV binding to the other NS5A conformations and even supported an asymmetrical Daclatasvir binding to the NS5A dimer.17 Our symmetrical and closed-conformation binding model is supported by the symmetrical nature of DCV and the symmetrical nature of the binding site we identified previously.16 In the current context, by a symmetrical ligand we strictly mean the symmetry of the pharmacophore of the compound, not the symmetry in their chemical composition. We show here that this symmetrical hypothesis fits a large range of DCV derivatives, which show great potency despite their different chemical structures. To further support this model, we describe here an in-depth computational analysis of the binding modes of six NS5A inhibitors. These inhibitors are shown in Figure 1 and include three symmetric (ABT-267;18

■

RESULTS AND DISCUSSION Mode of Binding of DCV and Analogues to HCV Genotype 1b NS5A. Different studies, including ours, have confirmed that NS5A inhibitors bind to the closed “back-toback” NS5A dimer.16,20 This is strongly supported by the symmetric features of DCV and many of its analogues (see Figure 1). The exact binding site of these inhibitors is located within close proximity to the dimer point of attachment to the membrane21 at a shallow solvent-exposed cleft that is distinctive to the NS5A dimer. This binding site is complemented from both sides by two N-terminal amphipathic protein helices that form the anchors to the ER membrane.21 It has been also suggested that the resistance patterns associated with a number of NS5A mutants result from disrupting the H-bond network of the surrounding residues; for example, between R30 and Q54 in GT-1b or between Q30 and R56 in GT-1a.22 However, given the exceptional picomolar potency of the currently available direct-acting NS5A inhibitors, most of these resistance mutations are not expected to result in ejection of the ligand from the binding site. An acceptable hypothesis is that these mutations enhance the flexibility of the NS5A dimer, enabling the attachment of NS5A dimer to its binding partners even in the presence of the bound ligands. Another hypothesis is that other protein conformations are populated in the mutant genotypes of the NS5A dimer and that current inhibitors do not target these conformations. For example, a recent crystal structure for the NS5A dimer of the GT-1a genotype (PDB ID 4CL1) 22 has shown that uncommon forms of NS5A dimerization are in fact possible (Figure 2). However, it is not clear if these forms are physiologically relevant. In our previously published model of DCV binding to NS5A, we showed that DCV binds in a symmetrical binding mode to the NS5A dimer.16 This symmetric binding mode is supported by hydrogen bonds from Arg30 and Gln54 from each monomer. This earlier model included a comprehensive

Figure 1. 2D chemical structures of the seven inhibitors under study. All potent direct-acting NS5A inhibitors are characterized by the presence of a central aromatic/hydrophobic biaryl core and two terminal peptide caps.

SYN-395 and SYN-7769) and three asymmetric (AZD;7 GSK;15 and Ledipasvir (LDV)19) molecular scaffolds. The most favorable binding mode of each inhibitor against the NS5A protein of the HCV GT-1a and GT-1b genotypes is identified and discussed in the light of available experimental and computational analyses data. The effect of known resistance-related mutations, including the L31V and Y93H mutations, is discussed. The current study enhances our understanding of the mode of action of this novel class of

Figure 2. Essential pharmacophoric features for potent directly acting NS5A inhibitors. These features are a central biaryl core (orange spheres), two H-bond acceptor sites (cyan spheres), and two hydrophobic sites (green spheres). The pharmacophoric sites are mapped over Daclatasvir (green sticks). 873



Article

Figure 3. 3D (upper panel) and 2D (lower panel) interaction diagrams of the most stable predicted binding modes of Daclatasvir to the (a) NS5A GT-1a (WT) and (b) NS5A GT-1b. Daclatasvir is stabilized by forming H-bonds with nearby residues, such as Gln30 and Arg56 in GT-1a and Arg30 and Gln54 in GT-1b. The central biaryl core extends over the surface at the middle of the shallow binding cleft characterizing the back-to-back dimer of NS5A.

the nearby lipophilic residues as well as a more robust fitting to the binding site. Figure 3 shows the predicted binding modes of DCV to NS5A encoded by HCV GT-1a and in GT-1b, where the surrounding pocket is shown as a surface representation. The 2D ligand interaction diagrams of DCV binding to the two genotypes are also given in the same figure. The binding modes of DCV to the two variants are very similar. The central biaryl core extends over the central region of the pocket, whereas the two side peptide caps engage with the nearby residues through H-bonds with Gln30 and Arg56 in GT-1a and with Arg30 and Gln54 in GT-1b. Similar interaction patterns for other DCV analogues have been observed in both genotypes. The predicted binding modes of the symmetric ligands (ABT, DCV, SYN-395, and SYN-776) to GT-1b are shown in Figure 4, whereas Figure 5 displays the binding modes of the asymmetric ligands (AZD, GSK, and LDV), also to GT-1b. Here the explicit H-bond interactions are shown in the 3D interaction diagram to illustrate the discrepancies between the different ligands at a higher resolution. As shown in Figure 4, all ligands bind almost symmetrically to the designated pocket where the biaryl cores extend horizontally over the binding cleft and the peptide caps engage with H-bonding to Arg30 and Gln54. The lack of the imidazole ring in ABT, which has the bulkiest central triphenylpyrrolidine core, seems to deprive the ligand of a stable

conformational analysis of the NS5A dimer in both its open and closed forms and incorporated many significant resistancerelated mutations described in the literature. The model was able to not only explain the effects of these mutations on the predicted activity of DCV but also elucidate the essential pharmacophoric features that are required for a potent NS5A inhibitor.16 Building upon that, we used the coordinates of our DCV− NS5A model to carry out pharmacophore-based docking simulations for the DCV-related structures shown in Figure 1. On the basis of these calculations, we found that without exception and despite these molecules having variant chemical compositions, all tested molecules satisfy the pharmacophore features of our earlier model (see below). Almost all NS5A inhibitors fall in a distinct region of chemical space and are composed of three main segments. The middle segment is a hydrophobic aromatic biaryl core, and the two terminal segments (peptide caps) are of variable sizes and nature and mostly contain a central amide group.23 The chemical structures of the compounds under study are shown in Figure 1. The bulkiness of the two terminal peptide caps seems to enhance the potency of the inhibitors and to retain some activity against relatively resistant mutants (see below).15 This enhanced potency provided by these bulky groups is presumably due to the formation of a better interaction with 874



Article

Figure 4. Detailed representation of the most stable predicted binding mode of the symmetric inhibitors under study, (a) ABT, (b) DCV, (c) SYN395, and (d) SYN-776, to the NS5A GT-1b (WT). The ligand is depicted in Licorice representation. The homodimer of the NS5A protein is colored silver. Hydrogen-bonding residues Arg30 and Gln54 are presented using a Licorice representation and colored (carbon is yellow, oxygen is red, nitrogen is blue, and hydrogen is white); ligand carbons are colored cyan. Hydrogen bonds are represented as dashed black lines.

Figure 5. Detailed representation of the most stable predicted binding mode of the nonsymmetric inhibitors under study, (a) AZD, (b) GSK, and (c) LDV, to the NS5A GT-1b (WT). The ligand is depicted in Licorice representation. The homodimer of the NS5A protein is colored silver. Hydrogen-bonding residues Arg30 and Gln54 are presented using a Licorice representation and colored (carbon is yellow, oxygen is red, nitrogen is blue, and hydrogen is white); ligand carbons are colored cyan. Hydrogen bonds are represented as dashed black lines. A small panel above GSK displays the burial of the 4-spiroketal ring.

875



Article

Figure 6. Most stable predicted binding mode of SYN-395 to the NS5A GT-1a and GT-1b (WT). As can be seen, SYN-395 interacts similarly with the two genotypes, forming H-bonds with Gln30 and Arg56 in GT-1a and with Arg30 and Gln54 in GT-1b.

molecular volumes of the binding site. This molecular volume reduction hinders the inhibitors from completely fitting into the designated binding site and also impairs the inhibitor’s ability to lock the NS5A protein in the closed conformation, which is not capable of binding to viral RNA.12 Thus, inhibitors with bulky substituents will compensate for this reduced molecular volume, causing a full occupation of the binding site. As shown in Figure 5b, the central core of GSK extends over the shallow binding cleft. The 4-spiroketal ring penetrates deeply into the small pocket formed by residues Ser25, Pro29, Ile52, Gly96, and Pro97. The binding mode is also supported by a strong charge-assisted H-bond with Lys26, Arg30, and Gln54. Figure 5b also shows a surface representation of the NS5A dimer complexed with GSK, where the 4-spiroketal group is deeply buried. This mode of binding slightly differs from a recently published model for GSK binding to the NS5A protein in which the 4-spiroketal group was facing upward.15 The fourth asymmetric inhibitor we studied is Ledipasvir (LDV). LDV is a pseudosymmetric molecule with a bulkier central aromatic hydrophobic core (difluorofluorene) compared to the DCV biphenyl core.25 Connecting the biphenyl core with a difluoromethylene bridge forms the central difluorofluorene ring. This enables LDV to form a robust interaction with residues in the shallow binding site. The most favorable binding mode of LDV to the NS5A GT-1b protein is shown in Figure 5c. The drug preserves all essential interactions with NS5A as in the case of DCV, including the H-bonds with Gln54 and Arg30 and the hydrophobic interaction with Leu27. LDV (EC50 = 4 pM) is almost 4 times more potent than DCV (EC50 = 15 pM).13 Substituting the biphenyl core in DCV by the difluorofluorene ring seems to enhance the overall rigidity of the molecules characterized by LDV.13 This enhanced rigidity as well as the incorporation of the azabicyclo ring systems in the peptide cap portion of the molecules seems to have a positive impact on the potency of the difluorofluorene core-containing inhibitors versus the biphenyl core-containing inhibitors, particularly against the NS5A GT-1b genotype.13 Comparing the Modes of Binding in GT-1b and GT1a. It was interesting to investigate the reasons behind the enhanced activity of many NS5A inhibitors to GT-1b as compared to GT-1a. Understanding these differences at the atomic level should lead to the discovery of inhibitors with more balanced activity against the two genotypes. We started

interaction with one of the Arg residues. In contrast to other inhibitors in the symmetric ligands category (see below), ABT exhibits a relatively tilted conformation within the binding site. Other inhibitors in this category, including DCV, SYN-395, and SYN-776, exhibit perfectly symmetrical binding modes, similar to what we have described in our previous study where Arg30 and Gln54 engages with strong H-bonds to the imidazole ring and the side peptide caps on both sides of the biphenyl core. SYN-395 and SYN-776 exhibit very similar binding modes against the NS5A GT-1a and NS5A GT-1b. Both drugs bear the essential pharmacophoric elements of DCV that we have identified previously to be essential for a properly functioning NS5A DAA inhibitor. SYN-395 has a relatively weak anti-HCV activity when used as a single treatment (EC50 = 214 nm against GT-1a wild type); however, the inhibitor was capable of restoring the impressive antiviral picomolar potency of DCV against the resistant mutants of different genotypes.9 For example, adding SYN-395 to DCV enhances the DCV activity against the GT-1a Y93N mutant from 339 to 0.13 nm, that is, a 2600-fold increase in potency in the HCV replicon assay. This cooperative effect translates to a substantial viral decline (>3 log10) in the infectious virus assay. A similar effect was observed when SYN-776 was combined with BMS-393.9 The binding modes of all asymmetric scaffolds are shown in Figure 5. As expected, these inhibitors exhibited slightly different interactions on both sides of the central core compared to symmetrical molecules. The smallest member in this category, AZD, does not fill the binding site completely, thus explaining the relatively weaker inhibitory potency against the NS5A genotype 1b (EC50 = 7000 pM).24 As shown in Figure 5a, the AZD amide bridge connects the two central biphenyl groups. This central amide group forms a strong Hbond with Thr95. The amide-cyclopropyl group forms Hbonds with Gln54, whereas the trifluoromethoxy group is Hbonded to Thr56. At one side of the dimer, Tyr93 forms a strong cation−pi stacking interaction with the protonated piperazine nitrogen. Incorporating a 4-spiro ketal motif to the main scaffold of DCV produces GSK (Figure 1).15 The most favorable binding mode of GSK against NS5A GT-1b is shown in Figure 5b. The design philosophy behind GSK was to retain the picomolar potency against the resistant mutants (L31V and Y93H).15 In general, L31V and Y93H mutations result in reducing the 876



Article

Q30, L31, P32, H54, and Y93 in NS5A GT-1a and L28, R30, L31, P32, F37, Q54, P58, and Y93 in GT-1b.26,28,29 In our previous study, we were able to map some of these resistance mutations over the DCV binding site and to correlate their effects to the observed activity of DCV in these mutations. The two primary resistance sites of mutations in NS5A HCV GT-1b are at residues L31 and Y93.11,29,30 As discussed in a very recent paper, it is estimated that >17.0% of all GT-1b infections are prone to developing resistance compared to the GT-1a infection (7.1%).31 This is presumably related to the fitness of the mutants and their ability to develop a stable infection as it has been shown in other studies that GT-1b mutants have a better replication level than those observed in GT-1a.28 For example, the replication level of the GT-1a L31M mutant was found to be approximately 55% of the wild-type replication level, whereas it was approximately 99% in the case of the GT1b L31M mutant. For these reasons, here we focused our mutational analysis mainly on the GT-1b variant. Interestingly, it has been observed that the in vitro replicon and the in vivo clinical data are strongly correlated,32 meaning that correlating the effect of mutations at the protein level to the clinical observations is indeed possible. Below, we try to understand the origin of these mutations and add more details of the mechanisms we presented in our earlier published model.16 For all inhibitors under study, the protein structures of the most energetically favorable binding pose in the GT-1b wild type were mutated to the GT-1b Y93H mutant and the GT-1b Y93H−L31V double mutant. Molecular dynamics simulations (3 ns) and MM/GBSA binding free energy calculations were carried out with the same parameter and settings as discussed earlier. For reasons that will become clear later, we have created a GT-1a Y93H model bound to DCV and SYN-395 and we performed MD simulation and binding energy calculations with the same parameters and settings. This gives a total of 16 complexes for the mutated models. The first step in our analysis was to quantify the mutationinduced reduction in the total MMGBSA binding energies of the inhibitors compared to the wild type, starting with the prototype inhibitor, DCV. It is well established from previous studies that DCV has a relatively lower barrier of resistance compared to other direct-acting NS5A inhibitors. For GT-1b, our results showed that indeed DCV exhibited one of the largest reductions in the total binding energies in the GT-1b Y93H−L31V double mutant (ΔΔGtotal = 21.5 kcal/mol). This significant reduction in the binding energy is predominated mainly by electrostatics (ΔΔGeel = 23.3 kcal/mol), supported by vdW effect (ΔΔGvdw = 13.2 kcal/mol) and partially compensated for by a reduced desolvation penalty (ΔΔGsolv = −15.01 kcal/mol). The reduction in the binding energy was less severe in the case of the GT-1b Y93H single mutant (11.8 kcal/mol). For the DCV binding to the GT-1a Y93H mutant, the reduction in binding energy was found to be 12.01 kcal/ mol. These data fit nicely to the reported resistance patterns for DCV that showed the fold resistance for DCV to increase in the following order: GT-1b Y93H (19-fold, EC50 increases from 2.6 to 49.2 pM) < GT-1a Y93H (5367-fold, EC50 increases from 6 to 32200 pM) < GT-1b Y93H−L31V (8336-fold, EC50 increases from 2.6 to 21674.5 pM).26,28 However, the small reduction in binding energy in the case of the Y93H mutants for both genotypes, which is given by 0.3 kcal/mol, is not sufficient to justify the huge difference in the fold resistance (5367-fold in GT-1a compared to 19-fold in GT-1b).

our analysis with DCV. As shown in Figure 3, the binding mode of DCV exhibited almost identical symmetric binding modes against both GT-1a and GT-1b. In in vitro replicon assays, the anti-HCV activity of DCV against GT-1b is approximately 5− 10 times more potent than the corresponding activity against GT-1a.26 A clear explanation for this is the relative strength of the H-bond formed between Arg30 in GT-1b and DCV as compared with the one formed with Gln30 in GT-1a. The MMPBSA binding energies for the most favorable binding modes of DCV against GT-1a and GT-1b were decomposed per residue. As expected, and given the inherently better electrostatic binding capabilities of Arg compared to Gln, the electrostatic contribution of Arg30 to the binding energy of DCV against NS5A GT-1b was approximately −17.6 kcal/mol. In GT-1a, the per-residue electrostatic contribution of the two Gln30 residues was estimated as −4.4 kcal/mol. The H-bond occupancies of Arg30 in GT-1b and of Gln30 in GT-1a were also analyzed and were found to be approximately 90% of the simulation time with an average value of 3.0 Å. The reduced electrostatic binding that results from the loss of Arg30 to Gln30 in Gt-1a (R30Q) seems to be partially compensated by a nearby Arg56 residue. Our per-residue binding energy decomposition analysis showed that Arg56 contributes favorably in the DCV binding to GT-1a (ΔΔGele ≈ −10.0 kcal/mol). In GT-1b, a Thr56 residue replaces Arg56. The context of this substitution will be discussed in more detail in the following sections. Analogously, the binding discrepancies of other inhibitors to the two genotypes under study (GT-1a and GT-1b) were examined. For example, Figure 6 displays the predicted binding modes of SYN-395 against GT-1a and GT-1b. As shown in Figure 6, the most favorable binding mode of SYN-395 engages with a very similar binding environment in the two genotypes, supporting our previous model for the binding of DCV as well as other direct-acting NS5A inhibitors. The central biphenyl core extends over the center of the shallow cleft over Thr95 residues. The peptide caps are H-bonded to Gln30, Arg56, and His54 in GT-1a and to Arg30, Gln54, and Lys26 in GT-1b. Other inhibitors under study form similar interactions with the two genotypes to various degrees. Effects of fhe Y93H and L31V Mutations. A major problem for this exceptionally potent class of direct-acting NS5A inhibitors is the possibility of developing resistance among treated patients.20,27,28 Particularly for DCV, several drug-induced mutations have been thoroughly described in the literature. A number of these resistance mutations confer several thousand-fold resistance to DCV treatment. For example, in GT-1b, a L31V−Y93H-linked mutation has been reported to cause an approximately 8336-fold resistance to DCV treatment.28 This represents a dramatic increase in the fold resistance that results from the single GT-1b L31V (23fold) or GT-1b Y93H (19-fold), implying a cooperative effect of the two mutations.28 Thus, investigating the structural implications of these mutations that led to this catastrophic pattern of resistance is of great interest for several industrial and academic research groups worldwide, including our group. This will greatly support the ongoing research efforts to develop potent inhibitors with high barriers to resistance. Certainly, the initial step for the rational design of new NS5A inhibitors that are active against resistant genotypes is to characterize these mutations at the atomic level. Almost all clinically relevant mutations are primarily located in the vicinity of the designated DCV binding site. These sites include M28, 877



Article

mutant, the reduction in binding energy compared to the wild type was −7.01 kcal/mol. To gain more structural insight regarding the aforementioned mutations, we have also analyzed the relaxed complexes after the MD simulation to see if any peculiar structural consequence of these mutations is taking place, particularly with the GT-1b genotype. In principle, these structural consequences are presumably similar for all inhibitors given the fact that we have carried out a pharmacophore-based docking in the beginning to get those complexes. In our previous study, it was clear that the L31V mutation induces a disruption in the nearby lipophilic residues, which enhances the flexibility of the amphipathic helices, discussed earlier. Mutation in the Y93 site, known to be the most prevalent in the GT-1b-resistant mutants, resulted in two structural consequences. The first consequence was that mutating Y93 with the shorter and smaller H93 reduces the vdW interaction with the lipophilic group of the bound inhibitor, for example, in DCV. In the current study, we have decomposed the binding energy on a per-residue basis. It was found that the Y93H mutation would result in a loss of approximately 2.1 kcal/mol in the total binding energy in the case of DCV binding to GT-1b Y93H mutant compared to the wild type. The second consequence was to deprive the protein from an important H-bond formed by the Y93 hydroxyl and a nearby proline residue (P29). In our current model, we found that there is also a possibility of some water molecules to bridge this interaction. As shown in Figure 8, two water molecules connect Y93 to P29. Such interaction is

Consistent with other experimental observations, this suggests that a rather complex structural consequence of the Y93H mutation is taking place in the case of the NS5A GT-1a genotype compared with the N55A GT-1b genotype.26,28 Another point to consider is the inherent flexibility of the two N-terminal amphipathic helices. In our structural models for DCV binding to the NS5A GT-1b wild type and mutants, we observed that the average RMSD over the MD simulation for GT-1b wild type was given as 2.94 Å for the full model and as 1.36 Å if we excluded the two amphipathic helices from the RMSD calculations (residues 1−25). However, for the GT-1b Y93H mutant, those values translated to 3.25 and 1.74 Å, respectively. For the GT-1b Y93H−L31V double mutant, the average RMSD for the full model was given as 3.27 and 2.1 Å without the helices. This means that the models exhibit a slightly higher flexibility in the amphipathic helix region for the mutant versus the NS5A wild type. This higher flexibility of the mutant in the amphipathic helix region could enable the NS5A protein to bind to the ER membrane (see above) even in the presence of the small molecule. Together with the reduction in binding energy, this enhanced flexibility could explain the observed resistance pattern, at least for DCV. Also for other inhibitors under study, we have compared the reduction in binding energies that result from the Y93H mutation and the Y93H−L31V double mutation. This reduction data is given as a bar chart in Figure 7. As shown

Figure 7. Bar chart representation for the reduction in MMGBSA binding energies for the seven inhibitors under study going from NS5A GT-1b wild-type to the GT-1b Y93H mutant (red bars) and from the NS5A GT-1b wild-type to the GT-1b Y93H−L31V double mutant. As shown, with the exception of AZD and SYN-395, the Y93H−L31V double mutant exhibited a larger effect on the binding energies of the inhibitors under study compared to the Y93H single mutant.

Figure 8. Bridging water molecules that H-bond Tyr93 and Pro29. Mutation of Tyr93 will affect the formation of this water bridge and result in disruption of this favorable network of water-mediated Hbonding. The homodimer of the NS5A protein is colored silver. Hydrogen bonds are represented as dashed black lines; water molecules are shown in vdW representation.

in the figure and with the exception of AZD and SYN-395, the binding energies of all inhibitors are more affected by the Y93H−L31V double mutation compared to the Y93H single mutation. GSK was the least affected by the Y93H single mutation, which can explain its ability to retain an excellent activity profile against this mutant in vitro and in vivo.15 DCV was the second most affected inhibitor by the Y93H−L31V double mutation and the fourth most affected inhibitor by the Y93H mutation. For SYN-395 binding to the GT-1a Y93H

of course absent in the Y93H mutant because histidine is shorter and cannot form this water-mediated interaction. As water is a very important parameter to drive the molecular recognition in biomolecules, we believe that this can add an additional piece of information in solving the puzzle of the observed resistance.33−35 The Puzzle of NS5A GT-1a Q30R. Another intriguing point that we tried to solve was that the NS5A GT-1a Q30R mutant seems to have a significant negative impact on the potency of DCV, causing approximately 1217-fold resistance 878



Article

(EC50 increased from 6 to 7300 pM).26,28 This slightly contradicts our assumption that a strong H-bond donor in this site should stabilize the interaction with the inhibitor, as in the case of GT-1b wild type where position 30 is actually an Arg residue. We observed that there exists a nearby Arg56 residue in GT-1a that is replaced by a neutral Thr56 residue in GT-1b. In GT-1a, Arg56 seems to interact favorably with DCV. According to our model, the presence of two strongly positively Arg residues (R30 and R56) in the vicinity of the DCV binding site might cause an electrostatic repulsion in the case of the Q30R mutant in GT-1a.

step required the use of the SHAKE algorithm43 to constrain all H-containing bonds. The temperature was fixed at 300 K using Langevin dynamics with the collision frequency γ of 1 ps−1. Pressure was maintained at 1 bar using the Berendsen barostat. Edge effect is removed by applying periodic boundary conditions. Long-range electrostatics are computed using the particle mesh Ewald (PME) summation method.44 All MD simulations were performed using the PMEMD module of AMBER14 using Compute Canada and IBM Blue Gene/Q supercomputer computational resources. Calculations of the Binding Free Energies. The binding free energies of the different inhibitors NS5A were calculated using the molecular mechanics/generalized Born surface area (MM/GBSA) module in AMBER over the 3 ns MD trajectory.45 Trajectory frames were saved every 2 ps spacing. For each protein−inhibitor complex, every second frame of the corresponding trajectory was used in the calculation, which creates a spacer of 4 ps between each frame. The free energy of each system was calculated as the sum of the gas phase molecular mechanics energy (EMM) and the solvation free energies (Gsolv):

■

METHODS Ligand Docking Simulations to the NS5A Protein. All ligands were prepared with low-energy tautomers and conformers generated using the LigPrep module as implemented in Schrodinger.36,37 For the protein structures, we used the atomic coordinates from our previously published model for DCV bound to NS5A GT-1b.16 Exhaustive docking simulations were carried out through the pharmacophore placement method in MOE.38,39 The previously developed pharmacophore for ligand binding to the GT-1b NS5A genotype was used. This pharmacophore, shown in Figure 2, consists of two central hydrophobic features, two hydrogen bond acceptors, represented by amide carbonyl and imidazole nitrogen in DCV, and a terminal bulky hydrophobic group, exemplified by a pyrollidine ring in DCV. All poses were rescored using the most recent GBVI/WSA dG38 scoring function as implemented in MOE. The 20 best-scoring poses for each docked ligand were retained after rescoring with the GBVI/WSA dG scoring function. For the GT-1a NS5A genotype, the atomic coordinates of the closed GT-1b NS5A genotype were converted to the GT-1a NS5A genotype by mutating binding site residues to their GT1a variants. We were also interested in analyzing the interaction of the ligands to two resistant mutants of the GT-1b NS5A genotype (the Y93H mutant and the Y93H−L31V double mutant).29,37 These mutants were obtained by mutating the corresponding protein residues of the most energetically favorable binding pose of the seven inhibitors under study to the GT-1b NS5A wild type. Molecular Dynamics Simulations. A maximum of 20 topscoring protein−ligand complexes for each ligand were selected for subsequent MD simulation in AMBER14.40 The AMBERff99SB force field for the protein and the general AMBER force field (GAFF) for the ligands were used.41 Ligand charges were calculated at the semiempirical AM1-BCC procedures in Antechamber.42 Each system was energy minimized for four 5000 minimization steps. During the first three minimization steps, protein and ligand atoms were constrained to their original positions by a force constant of 100, 50, 5 kcal mol−1 Å−2, respectively, and then freely minimized for an additional 5000 steps. Each system was heated to 300 K at the NVT ensemble for 50 ps with protein, and ligand heavy atoms were constrained to their original positions by a force constant of 5 kcal mol−1 Å−2 and using a 1 fs time step to integrate the equations of motion. The systems were then equilibrated at the NPT ensemble for four equilibration steps of 10 ps each with gradual relief of the constraints from protein backbone atoms and finally equilibrated for 25 ps free MD simulation using a 2 fs time step. This was followed by 6 × 0.5 ns production MD simulation using a 2 fs integration time step, giving a total of 3 ns simulation trajectories. The use of a 2 fs integration time

G = EMM + Gsolv

The MM energies were calculated using the SANDER program in AMBER, and the solvation free energy part was calculated using the generalized-born (GB) solvation model to approximate the explicit solvent. Binding free energies were estimated as the average differences of the free energies between the protein (P) and the ligand (L) in their bound (PL) and unbound states:46 ΔGbinding = G PL − G P − G L

Hydrogen bonding occupancies were extracted from the 3 ns production simulation trajectory using the CPPTRAJ program in AMBER. A distance of 3.5 Å and an angle of 120° were used as cutoff values for accepting the hydrogen bonds.

■

CONCLUSION Hepatitis C virus is still a major and global health problem. The current standard of care involve drugs that directly target the viral proteins, known as direct-acting antivirals. NS5A represents a very promising HCV target, and drugs binding to this protein exhibit exceptional antiviral efficacy. Despite this exceptional potency, very little is known about how NS5A inhibitors interact with the protein. Daclatasvir is the most studied DAA targeting NS5A. Our earlier study identified the binding site of DCV in NS5A and revealed a symmetrical mode of binding of the drug. Here, we have built upon this model and investigated the modes of binding of a number of DCV-like structures in two different genotypes, namely, GT-1b and GT-1a. The NS5A inhibitors we studied included both symmetric and asymmetric molecular scaffolds. Both types of molecules satisfied our previously described pharmacophore, although the two classes revealed a slightly different mode of binding. Additionally, it was clear that bulky derivatives achieve a more robust interaction through engagement with small side pockets in the proximity of the main shallow pocket at the NS5A dimer interface. Our binding energy calculations correlated very well with the observed in vitro replicon and mutational data. Substitution of Tyr93 residues with other residues in the mutant genotypes resulted in inhibiting the formation of H879



Article

(12) Ascher, D. B., Wielens, J., Nero, T. L., Doughty, L., Morton, C. J., and Parker, M. W. (2014) Potent hepatitis C inhibitors bind directly to NS5A and reduce its affinity for RNA. Sci. Rep. 4, 4765. (13) Kwon, H. J., Xing, W., Chan, K., Niedziela-Majka, A., Brendza, K. M., Kirschberg, T., Kato, D., Link, J. O., Cheng, G., Liu, X., and Sakowicz, R. (2015) Direct binding of ledipasvir to HCV NS5A: mechanism of resistance to an HCV antiviral agent. PLoS One 10, e0122844. (14) Zillner, K., Jerabek-Willemsen, M., Duhr, S., Braun, D., Langst, G., and Baaske, P. (2012) Microscale thermophoresis as a sensitive method to quantify protein: nucleic acid interactions in solution. Methods Mol. Biol. 815, 241−252. (15) Kazmierski, W. M., Maynard, A., Duan, M., Baskaran, S., Botyanszki, J., Crosby, R., Dickerson, S., Tallant, M., Grimes, R., Hamatake, R., Leivers, M., Roberts, C. D., and Walker, J. (2014) Novel spiroketal pyrrolidine GSK2336805 potently inhibits key hepatitis C virus genotype 1b mutants: from lead to clinical compound. J. Med. Chem. 57, 2058−2073. (16) Barakat, K. H., Anwar-Mohamed, A., Tuszynski, J. A., Robins, M. J., Tyrrell, D. L., and Houghton, M. (2015) A Refined Model of the HCV NS5A protein bound to daclatasvir explains drug-resistant mutations and activity against divergent genotypes. J. Chem. Inf. Model. 55, 362−373. (17) Nettles, J. H., Stanton, R. A., Broyde, J., Amblard, F., Zhang, H., Zhou, L., Shi, J., McBrayer, T. R., Whitaker, T., Coats, S. J., Kohler, J. J., and Schinazi, R. F. (2014) Asymmetric binding to NS5A by daclatasvir (BMS-790052) and analogs suggests two novel modes of HCV inhibition. J. Med. Chem. 57, 10031−10043. (18) DeGoey, D. A., Randolph, J. T., Liu, D., Pratt, J., Hutchins, C., Donner, P., Krueger, A. C., Matulenko, M., Patel, S., Motter, C. E., Nelson, L., Keddy, R., Tufano, M., Caspi, D. D., Krishnan, P., Mistry, N., Koev, G., Reisch, T. J., Mondal, R., Pilot-Matias, T., Gao, Y., Beno, D. W., Maring, C. J., Molla, A., Dumas, E., Campbell, A., Williams, L., Collins, C., Wagner, R., and Kati, W. M. (2014) Discovery of ABT267, a pan-genotypic inhibitor of HCV NS5A. J. Med. Chem. 57, 2047−2057. (19) Link, J. O., Taylor, J. G., Xu, L., Mitchell, M., Guo, H., Liu, H., Kato, D., Kirschberg, T., Sun, J., Squires, N., Parrish, J., Keller, T., Yang, Z. Y., Yang, C., Matles, M., Wang, Y., Wang, K., Cheng, G., Tian, Y., Mogalian, E., Mondou, E., Cornpropst, M., Perry, J., and Desai, M. C. (2014) Discovery of ledipasvir (GS-5885): a potent, once-daily oral NS5A inhibitor for the treatment of hepatitis C virus infection. J. Med. Chem. 57, 2033−2046. (20) Issur, M., and Gotte, M. (2014) Resistance patterns associated with HCV NS5A inhibitors provide limited insight into drug binding. Viruses 6, 4227−4241. (21) Penin, F., Brass, V., Appel, N., Ramboarina, S., Montserret, R., Ficheux, D., Blum, H. E., Bartenschlager, R., and Moradpour, D. (2004) Structure and function of the membrane anchor domain of hepatitis C virus nonstructural protein 5A. J. Biol. Chem. 279, 40835− 40843. (22) Lambert, S. M., Langley, D. R., Garnett, J. A., Angell, R., Hedgethorne, K., Meanwell, N. A., and Matthews, S. J. (2014) The crystal structure of NS5A domain 1 from genotype 1a reveals new clues to the mechanism of action for dimeric HCV inhibitors. Protein Sci. 23, 723−734. (23) Cordek, D. G., Bechtel, J. T., Maynard, A. T., Kazmierski, W. M., and Cameron, C. E. (2011) Targeting the Ns5a protein of Hcv: an emerging option. Drugs Future 36, 691−711. (24) Gish, R. G., and Meanwell, N. A. (2011) The NS5A replication complex inhibitors: difference makers? Clin. Liver. Dis. 15, 627−639. (25) Lawitz, E. J., Gruener, D., Hill, J. M., Marbury, T., Moorehead, L., Mathias, A., Cheng, G., Link, J. O., Wong, K. A., Mo, H., McHutchison, J. G., and Brainard, D. M. (2012) A phase 1, randomized, placebo-controlled, 3-day, dose-ranging study of GS5885, an NS5A inhibitor, in patients with genotype 1 hepatitis C. J. Hepatol. 57, 24−31. (26) Gao, M., Nettles, R. E., Belema, M., Snyder, L. B., Nguyen, V. N., Fridell, R. A., Serrano-Wu, M. H., Langley, D. R., Sun, J. H.,

bonds with water molecules and a decrease in the free energy of binding to ligands. We believe that the data presented in this paper will facilitate ongoing efforts toward enhancing the efficacy of DCV-like drugs.

■

AUTHOR INFORMATION

Corresponding Author

*(K.B.) E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Li Ka Shing Institute of Virology and the Li Ka Shing Applied Virology Institute at the University of Alberta, the Alberta Cancer Foundation, and The Natural Sciences and Engineering Research Council of Canada (NSERC). All simulations were made possible through the Compute Canada and IBM Blue Gene/Q supercomputer computational resources.

■

REFERENCES

(1) Pawlotsky, J. M. (2013) Hepatitis C virus: standard-of-care treatment. Adv. Pharmacol. 67, 169−215. (2) Gane, E. (2012) Future perspectives: towards interferon-free regimens for HCV. Antiviral Ther. 17, 1201−1210. (3) Deuffic-Burban, S., Schwarzinger, M., Obach, D., Mallet, V., Pol, S., Pageaux, G. P., Canva, V., Deltenre, P., Roudot-Thoraval, F., Larrey, D., Dhumeaux, D., Mathurin, P., and Yazdanpanah, Y. (2014) Should we await IFN-free regimens to treat HCV genotype 1 treatment-naive patients? A cost-effectiveness analysis (ANRS 95141). J. Hepatol. 61, 7−14. (4) Conteduca, V., Sansonno, D., Russi, S., Pavone, F., and Dammacco, F. (2014) Therapy of chronic hepatitis C virus infection in the era of direct-acting and host-targeting antiviral agents. J. Infect. 68, 1−20. (5) Asselah, T., and Marcellin, P. (2011) New direct-acting antivirals’ combination for the treatment of chronic hepatitis C. Liver Int. 31 (Suppl. 1), 68−77. (6) O’Boyle, D. R., 2nd, Nower, P. T., Gao, M., Fridell, R., Wang, C., Hewawasam, P., Lopez, O., Tu, Y., Meanwell, N. A., Belema, M., Roberts, S. B., Cockett, M., and Sun, J. H. (2016) Synergistic activity of combined NS5A inhibitors. Antimicrob. Agents Chemother. 60, 1573−1583. (7) Kohler, J. J., Nettles, J. H., Amblard, F., Hurwitz, S. J., Bassit, L., Stanton, R. A., Ehteshami, M., and Schinazi, R. F. (2014) Approaches to hepatitis C treatment and cure using NS5A inhibitors. Infect. Drug Resist. 7, 41−56. (8) Wang, C., Jia, L., O’Boyle, D. R., 2nd, Sun, J. H., Rigat, K., Valera, L., Nower, P., Huang, X., Kienzle, B., Roberts, S., Gao, M., and Fridell, R. A. (2014) Comparison of daclatasvir resistance barriers on NS5A from hepatitis C virus genotypes 1 to 6: implications for crossgenotype activity. Antimicrob. Agents Chemother. 58, 5155−5163. (9) Sun, J. H., O’Boyle, D. R., 2nd, Fridell, R. A., Langley, D. R., Wang, C., Roberts, S. B., Nower, P., Johnson, B. M., Moulin, F., Nophsker, M. J., Wang, Y. K., Liu, M., Rigat, K., Tu, Y., Hewawasam, P., Kadow, J., Meanwell, N. A., Cockett, M., Lemm, J. A., Kramer, M., Belema, M., and Gao, M. (2015) Resensitizing daclatasvir-resistant hepatitis C variants by allosteric modulation of NS5A. Nature 527, 245−248. (10) Bartenschlager, R., Lohmann, V., and Penin, F. (2013) The molecular and structural basis of advanced antiviral therapy for hepatitis C virus infection. Nat. Rev. Microbiol. 11, 482−496. (11) Nakamoto, S., Kanda, T., Wu, S., Shirasawa, H., and Yokosuka, O. (2014) Hepatitis C virus NS5A inhibitors and drug resistance mutations. World J. Gastroenterol. 20, 2902−2912. 880



Article

O’Boyle, D. R., 2nd, Lemm, J. A., Wang, C., Knipe, J. O., Chien, C., Colonno, R. J., Grasela, D. M., Meanwell, N. A., and Hamann, L. G. (2010) Chemical genetics strategy identifies an HCV NS5A inhibitor with a potent clinical effect. Nature 465, 96−100. (27) Poveda, E., Wyles, D. L., Mena, A., Pedreira, J. D., CastroIglesias, A., and Cachay, E. (2014) Update on hepatitis C virus resistance to direct-acting antiviral agents. Antiviral Res. 108, 181−191. (28) Fridell, R. A., Qiu, D., Wang, C., Valera, L., and Gao, M. (2010) Resistance analysis of the hepatitis C virus NS5A inhibitor BMS790052 in an in vitro replicon system. Antimicrob. Agents Chemother. 54, 3641−3650. (29) Yoshimi, S., Ochi, H., Murakami, E., Uchida, T., Kan, H., Akamatsu, S., Hayes, C. N., Abe, H., Miki, D., Hiraga, N., Imamura, M., Aikata, H., and Chayama, K. (2015) Rapid, sensitive, and accurate evaluation of drug resistant mutant (NS5A-Y93H) strain frequency in genotype 1b HCV by invader assay. PLoS One 10, e0130022. (30) Wang, C., Huang, H., Valera, L., Sun, J. H., O’Boyle, D. R., 2nd, Nower, P. T., Jia, L., Qiu, D., Huang, X., Altaf, A., Gao, M., and Fridell, R. A. (2012) Hepatitis C virus RNA elimination and development of resistance in replicon cells treated with BMS-790052. Antimicrob. Agents Chemother. 56, 1350−1358. (31) Dietz, J., Susser, S., Berkowski, C., Perner, D., Zeuzem, S., and Sarrazin, C. (2015) Consideration of viral resistance for optimization of direct antiviral therapy of hepatitis C virus genotype 1-infected patients. PLoS One 10, e0134395. (32) Gao, M. (2013) Antiviral activity and resistance of HCV NS5A replication complex inhibitors. Curr. Opin. Virol. 3, 514−520. (33) Robinson, C. R., and Sligar, S. G. (1993) Molecular recognition mediated by bound water. A mechanism for star activity of the restriction endonuclease EcoRI. J. Mol. Biol. 234, 302−306. (34) Levy, Y., and Onuchic, J. N. (2006) Water mediation in protein folding and molecular recognition. Annu. Rev. Biophys. Biomol. Struct. 35, 389−415. (35) Li, M., Hoeck, C., Schoffelen, S., Gotfredsen, C. H., and Meldal, M. (2016) Specific electrostatic molecular recognition in water. Chem.−Eur. J. 22, 7206−7214. (36) Sastry, G. M., Adzhigirey, M., Day, T., Annabhimoju, R., and Sherman, W. (2013) Protein and ligand preparation: parameters, protocols, and influence on virtual screening enrichments. J. Comput.Aided Mol. Des. 27, 221−234. (37) Wyles, D. L. (2013) Antiviral resistance and the future landscape of hepatitis C virus infection therapy. J. Infect. Dis. 207 (Suppl. 1), S33−S39. (38) Corbeil, C. R., Williams, C. I., and Labute, P. (2012) Variability in docking success rates due to dataset preparation. J. Comput.-Aided Mol. Des. 26, 775−786. (39) Labute, P. (2005) On the perception of molecules from 3D atomic coordinates. J. Chem. Inf. Model. 45, 215−221. (40) Case, D. A., Darden, T. A., Cheatham, T. E., Simmerling, C. L., Wang, J., Duke, R. E., Luo, R., Walker, R. C., Zhang, W., Merz, K. M., Roberts, B., Hayik, S., Roitberg, A., Seabra, G., Swails, J., Goetz, A. W., Kolossváry, I., Wong, K. F., Paesani, F., Vanicek, J., Wolf, R. M., Liu, J., Wu, X., Brozell, S. R., Steinbrecher, T., Gohlke, H., Cai, Q., Ye, X., Wang, J., Hsieh, M. J., Cui, G., Roe, D. R., Mathews, D. H., Seetin, M. G., Salomon-Ferrer, R., Sagui, C., Babin, V., Luchko, T., Gusarov, S., Kovalenko, A., and Kollman, P. A. (2014) AMBER 14, University of California, San Francisco, CA, USA. (41) Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004) Development and testing of a general amber force field. J. Comput. Chem. 25, 1157−1174. (42) Wang, J., Wang, W., Kollman, P. A., and Case, D. A. (2006) Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graphics Modell. 25, 247−260. (43) Ryckaert, J.-P., Ciccotti, G., and Berendsen, H. J. C. (1977) Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327−341.

(44) Darden, T., York, D., and Pedersen, L. (1993) Particle mesh Ewald: An N.log(N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089−10092 citeulike-article-id:297086.. (45) Miller, B. R., McGee, T. D., Swails, J. M., Homeyer, N., Gohlke, H., and Roitberg, A. E. (2012) MMPBSA.py: an efficient program for end-state free energy calculations. J. Chem. Theory Comput. 8, 3314− 3321. (46) Genheden, S., and Ryde, U. (2015) The MM/PBSA and MM/ GBSA methods to estimate ligand-binding affinities. Expert Opin. Drug Discovery 10, 449−461.

881


A Comprehensive Computational Analysis for the Binding Modes of

Recommend Documents