Asymmetric binding to NS5A by daclatasvir (BMS-790052) and analogs suggests two novel modes of HCV inhibition.

Symmetric, dimeric daclatasvir (BMS-790052) is the clinical lead for a class of picomolar inhibitors of HCV replication. While specific, resistance-bearing mutations at positions 31 and 93 of domain I strongly suggest the viral NS5A as target, structural mechanism(s) for the drugs' activities and resistance remains unclear. Several previous models suggested symmetric binding modes relative to the homodimeric target; however, none can fully explain SAR details for this class. We present semiautomated workflows to model potential receptor conformations for docking. Surprisingly, ranking docked hits with our library-derived 3D-pharmacophore revealed two distinct asymmetric binding modes, at a conserved poly-proline region between 31 and 93, consistent with SAR. Interfering with protein-protein interactions at this membrane interface can explain potent inhibition of replication-complex formation, resistance, effects on lipid droplet distribution, and virion release. These detailed interaction models and proposed mechanisms of action will allow structure-based design of new NS5A directed compounds with higher barriers to HCV resistance.

Introduction

Hepatitis C virus (HCV) infection is a global epidemic with associated high risk for serious liver disease.[1] Compound 1 (daclatasvir, BMS-790052) is the leading representative of a new class of direct-acting antiviral agents (DAA) against HCV infection that target the viral nonstructural protein 5A (NS5A). This family of compounds includes some of the most active antiviral compounds tested, with low picomolar median effective concentration (EC50) in HCV replicon assays.[2−5] Three structurally related compounds currently in clinical trials, 1, 2 (GSK-2336805), and 3 (GS-5885), are illustrated in Chart 1. Because NS5A lacks known enzymatic activity, the specific mechanism(s) for the extraordinary potency of this class of antiviral drugs is not yet clear. While cell-based studies have shown that NS5A is critical for viral replication,[6−8] clinical studies suggest these drugs inhibit multiple stages of viral release.[9,10] Most recently, NS5A-DAA have been shown to directly disrupt formation of the membranous viral replication complexes.[11]

Chart 1

Structurally Similar NS5A Directed Inhibitors Currently in Clinical Trialsa

The compounds 1 (BMS-790052), 2 (GSK-2336805), 3 (GS-5885) share two peptidic caps connected via an aromatic linker and are thought to bind the same site on the NS5A protein.

All reported NS5A-DAA rapidly select for multiple genotype-specific mutations in NS5A that markedly reduce efficacy. For example, in genotype 1b (Gt1b), a single mutation of L31 V or Y93H imparts 28- or 24-fold resistance to 1, respectively. However, the double mutation (31/93) imparts over 14 000-fold resistance in vitro (Table 1).[4] In clinical trials, compound 1 caused a rapid drop in viremia in responders but selected for the same 31/93 mutations in subjects with persistent Gt1b-infections.[2,12,13]

Table 1

In Vitro Genotype 1b Replicon Activity/Resistance Profile of Daclatasvir 1 Used for Structural Modeling Designa

substitution, genotype 1b	replication level (%), average ± SD	EC50 (pM), average ± SD	fold resistance
WT	100	2.6 ± 0.9	1
L31M	99 ± 23	8.4 ± 1.9	3
L31 V	158 ± 54	71.7 ± 20.3	28
Y93H	27 ± 16	62.3 ± 24.4	24
L31M + Y93H	70 ± 68	18271 ± 16512	7105
L31V + Y93H	50 ± 38	38031 ± 33429	14789

Bolded fold resistance values highlight the effect of double position vs single mutations at positions 31 and 93 on drug sensitivity, suggesting a shared role in drug binding. Data adapted with permission from (ref (4)).

Recently reported medicinal chemistry efforts have optimized early high throughput screening hits into potent clinical leads with improved pan-genotypic activity and resistance profiles.[14−17] However, further structure-based analysis/optimization has been hindered by the lack of crystallographic complexes with bound drug or theoretical binding-site models that explained the general structure–activity relationships (SARs) across the various series of NS5A compounds. The goal of this work was to develop receptor models for evaluating structural features of our NS5A-DAA analog series associated with either high activity change or sensitivity to resistance from mutations in NS5A domain I (D-I) near residues 31 and 93. Our approach was to systematically explore multiple theoretical binding models that could be filtered and validated using experimental activities from our NS5A-DAA analog library.

Current Structural Models of NS5A-DAA Binding

The 3D structure of NS5A associated with DAA binding and activity is controversial with six groups presenting different interaction models.[15,18−20] Groups at Bristol-Myers Squibb (BMS) and GlaxoSmithKline (GSK) each presented a unique theoretical model of the DAA binding to D-I based upon different experimental homodimers of truncated Gt1b-NS5A to map their drug’s interactions.[15,20] A recent theoretical study used large-scale molecular dynamics (MD) simulations to sample local conformational space based on the same two X-ray derived dimer forms followed by automated docking across the multiple receptor states to locate an energy “refined” binding mode.[19] The fourth group compared experimental NS5A RNA binding activity to theoretical docking against modeled Gt1b dimers.[21] A fifth model was derived through pharmacophore-based docking to a new homodimeric form observed in crystal structures of genotype 1a (Gt1a) NS5A.[18] All five groups suggested that the NS5A-DAAs bind symmetrically near the core homodimer interface including Y93 of each monomer, but the site interaction details differ greatly depending upon the dimeric forms of the receptor and methods used. Unlike these previous works, the most recent study finds primarily asymmetric binding modes when docking was preformed without theoretical alignments of the membrane-binding region that is missing from the experimental Gt1b dimer forms.[11]

In line with these discrepancies, a recent review of HCV proteins describes NS5A as “one of the most enigmatic” drug targets.[22] Because NS5A has no known enzymatic activity, there is no specific active site region to model. Structurally, the NS5A protein is 447 residues in length and divided into three domains (Figure 1a). D-I and portions of domain II appear to be required for viral replication, while domain III may be required for viral packaging and release.[23,24] D-I is the only domain with inherent structure and is characterized by its membrane-binding and cytosolic subdomains.[25,26] Sequence alignments of three solved D-I structural fragments that highlight positions of drug resistance residues 31 and 93 are shown in Figure 1b (PDB code 1R7G,[9] residues 1–31; PDB code 1ZH1,[25] residues 36–198; PDB code 3FQQ,[27] residues 32–191).

Figure 1

Three-domain structure of NS5A and alignment of domain 1 (D-I) sequences used for experimental structure evaluation and 3D modeling. (a) Three domain structure of full-length NS5A. (b) Domain 1 (D-I): subdomain detail and sequence alignments of HCV genotypes 1a and 1b with experimental structures from NMR (1R7G) and X-ray crystallography (1ZH1, 3FQQ). Primary resistance-bearing mutation positions 31 and 93 from Table 1 are boxed and highlighted in red. AH-D-Ia linker containing a conserved PxxPxxP motif is boxed orange. There are no experimental coordinates for a complete AH-linker.

The N-terminal region of D-I is required for biological membrane-binding of NS5A; however, that function also has prevented experimental structural characterization of full-length D-I in vitro. The isolated N-terminal 31 residues of genotype 1a (Gt1a) were characterized by elegant protein NMR (PDB code 1R7G)[9] and biochemical studies as a flexible amphipathic helical region (AH) that affects intracellular localization of the protein and is critical for HCV replication.[7,28,29] The NMR models of AH position the α-helix “in-plane” with the membrane surface having polar groups toward the cytosol and hydrophobic groups toward the lumen.[7]

The cytosolic portion of Gt1b D-I was solved with X-ray diffraction by two independent laboratories, but each used different truncated sequences and produced alternative homodimeric 3D models.[25,27] Tellinhausen et al. used fragments with peptide residues 25–198 that revealed a novel, ordered Zn2+ binding domain (D-Ia) between residues 36–100 that coordinated a zinc atom via four cysteines.[25] They also found a putative RNA binding subdomain 1b (D-Ib) in region 101–213 that is formed at the homodimer interface (PDB code 1ZH1);[25] however, positions of residues 25–35 that would link cytosolic D-Ia to membrane-bound AH could not be determined.[25] Love et al. used shorter peptide fragments including residues 32–191 that resulted in monomeric 3D structures of D-Ia and D-Ib subdomains that superimpose with the Tellinhausen et al. structure but pack in a different dimeric form. While the Love et al. structure had an identical D-Ia Zn2+ binding region, the D-Ib subdomains were not aligned for RNA binding because of the different dimer configuration (PDB code 3FQQ).[27] The two dimeric forms have complicated the building and interpretation of full-length D-I models that link AH to D-Ia/D-Ib.[25,27] For example, Bartenschlager et al. showed that the AH from NMR (containing residue L31) can be manually aligned close to Y93 residues at the interface of D-Ia/D-Ib X-ray derived dimer form by orienting the N-terminal end outward; however, this creates potential site models with very different characteristics.[22] Similar to Bartenschlager et al., all other published models of NS5A-DAA align the AH with the N-terminus extended away from each symmetric dimeric core form; however, such methods of alignment have been described as “arbitrary”.[11]

Herein, an alternative approach to obtain D-I models that systematically explores a broad range of AH-D-Ia/D-Ib receptor conformations available to each of the dimer forms and a method to evaluate their biological relevance is presented. We discovered a low-energy conformation for a conserved poly-proline motif (PxxPxxP) in region 29–35 with the potential to position in-plane with a membrane surface, consistent with its biological function linking AH with D-Ia/D-Ib. This low-energy PxxPxxP linker conformation also minimizes distances between highly mutable residues L31 and Y93. Proline rich motifs, such as PxxPxxP on the NS5A linker region, are often involved in protein interactions and signaling, so we hypothesized that this region could be involved in drug binding.[30] By using experimentally solved membrane-binding proteins with homologous PxxPxxP regions as “flexible” AH modeling templates, we built a family of full-length Gt1b NS5A-D-I receptor models that revealed a novel packing of AH to the D-Ia/D-Ib dimer fold that was low in chemical strain, minimized distances between residues L31 and Y93, lies in plane with the membrane surface, and was suitable for automated docking of NS5A-DAA across a range of receptor conformations. All models were evaluated using empirical data from our library of 64 NS5A-DAA analogs for pharmacophore elucidation, ranking of docked results, and subsequent testing of SAR within the multiple receptor conformations. The results suggest two novel asymmetric drug binding modes that provide structural explanation for selected resistance and SAR and offer a two-step mechanism of action to explain the broad activity profiles observed across compounds of this class.

Results and Discussion

Our NS5A-directed drug discovery program included parallel computational modeling and medicinal chemistry approaches to develop tools for elucidation of mechanism(s) explaining the picomolar activity of 1 and other compounds in this class. Specifically, our computational methods utilized publicly available structural data and focused on regions associated with high drug resistance for development of potential binding models. Our chemistry group designed analogs that were systematically evaluated for structure–activity relations to help optimize our models and develop improved compounds toward clinical advancement. The chemistry efforts, anti-HCV activity, and SAR have been discussed previously,[31−33] and this work focuses on development of chemical and biological modeling methods and their potential use in understanding the patterns revealed by drug activities observed in our Gt1b replicon assay.

Structure based sequence alignments are presented in Figure 1b. High resistance for combined mutations at positions 31 and 93 for Gt1b suggests a shared role in drug binding (Table 1), and the residues are highlighted in red boxes. Unfortunately, no single experimental study provided contiguous structural data for the D-I linker region (residues 26–36) that contains residue 31 and connects the cytosolic Zn2+ binding subdomain to the membrane-bound AH because of disorder in the crystal for that region of the protein.

3D Alignments and Experimental Structure Analysis

Three-dimensional coordinates were downloaded from the RCSB Protein Data Bank (PDB) and compared as described in Experimental Methods. The Gt1b D-I cytosolic fragments from the previously described X-ray diffraction studies were evaluated for sequence and structural variances between the two homodimers forms. PDB code 1ZH1, from Tellinhausen et al., included coordinates for residues 36–198;[25]3FQM and 3FQQ(27) structures, by Love et al., include residues 32–191.[27] It is important to note that the protein for the Tellinhausen et al. structure (solved in the laboratory of Charles Rice), 1ZH1, included the residues 25–35 of the flexible linker region, but they were disordered in the crystal and positions could not be incorporated in the model.[25] While the 3D monomeric subunits aligned with a root mean squared deviation (rmsd) of only 0.5 Å between 154 concurrent α-carbons of the peptide backbone, the position of the dimeric subunits varied greatly relative to one another (Figure S-1, Ia). The most divergent residues between the aligned monomers were between F36 and F37 with α-carbon distances of 6.2 and 9.2 Å, respectively, and may play a role in the different dimeric packing. Since the physiological relevance of neither dimeric conformation has been determined, we tested each as a template for manual and automated linking of the membrane binding AH subdomains. Our manual fitting results were similar to recently published models;[15,20] however, no energy-minimized model remained consistent with planar AH/membrane alignment as deduced by NMR. To assess a potential role for geometry of the conserved PxxPxxP motif related to the linking of the D-Ia core to a membrane-bound AH, we defined eight residues 28–36 as our flexible “linker” and used Modeller[34] to generate 20 diverse, low-energy peptide conformations for each dimer form (10 per monomer) (Figure S-1, Ib). We then ranked each conformation for minimum distances between C-alpha of residues 31 and 93 from either monomer as described in methods (Figure S-1, Ic). Notably, the lowest 31–93 inter-residue distance (9.4 Å) of the linker models that had planar alignment was found between residues located on different subunits of the 1ZH1 form. Appending this novel peptide linker model to each 1ZH1 monomer extended the dimer interface from position Y93 to P97 and resulted in a potential binding region having symmetry and shape/size similarity to compound 1 suitable for further model building (Figure S-1, Id). It was interesting to note that residues Y93 to P97 also formed the alternative interface for the 3FQM and 3FQQ dimers, suggesting a valid protein–protein interaction surface.[27]

Fitting Homologous Region of Sorting Nexin 5 (SNX5) Phox (PX) Domain to NS5A PxxPxxP Linker Template

To evaluate further biological relevance our theoretical linker conformation, we searched for proteins having experimentally solved PxxPxxP motifs and identified two structures of the sorting nexin 5 (SNX5) phox (PX) domain (PDB codes 3HPB and 3HPC).[35] SNX5 has been shown to influence membrane curvature and receptor internalization similar to NS5A’s attributed functions.[36] Both of the experimental SNX5 structures fit the theoretical NS5A PxxPxxP linker with a root-mean-square deviation of <1.5 Å for the prolines. Surprisingly, each experimental SNX5 structure contained a 180° turn (adjacent to the PxxPxxP motif) that aligns the N-terminal amphiphilic helices toward the NS5A dimer core, which is unlike any previously reported model of NS5A-DI (Figure S-1, IIb). Using the trimmed SNX5 region to align the NMR models of NS5A residues 1–31 resulted in a fit of backbone atoms with rmsd < 1 Å that also positioned key side chain residues, including the four conserved tryptophans (4, 9, 11, and 22), consistent with membrane alignments observed by NMR and reported by Penin et al.[7] (Figure S-1, IIc,d). We next tested this hybrid NS5A/SNX5 structure as a template for sequence-based homology modeling of the full-length NS5A domain 1.

Homology Modeling of Full-Length NS5A Domain 1

We generated families of models for AH connected NS5A domain 1 using the SNX5-NS5A chimera with fixed and flexible protocols described in Experimental Methods that resulted in 10 low-energy, structurally diverse potential receptors for each set (Figure S-1, III).

SAR Library Development

A Pipeline Pilot workflow was used to format our 68 previously reported NS5A analogs with Gt1b replicon EC50 values ranging from 2 pM to >106 pM for SAR and pharmacophore analysis in MOE (Figure S-2). The compound library systematically explores structural variations to three scaffold families that share a proline–valine–carbamate cap[31−33] (Figure 2).

Figure 2

NS5A inhibitors: series I, II, III and common pharmacophore elucidation. Best model suggests a specific monomeric conformation with two aromatic features (gold spheres) and two buried donor/receptor pairs (blue spheres) aligned relative to a peptide-like “γ-turn” (blue arrow) that is shared by all picomolar actives across the three chemical series. Steepest activity cliffs were selected from each series as described in methods section and Figure S-2.

Pharmacophore Elucidation

A subset of 23 highly active compounds (picomolar EC50) was identified across the three chemical series for pharmacophore elucidation as described in Experimental Methods. Surprisingly, our best model suggested the primary pharmacophore responsible for high activity of the symmetric, dimeric compound 1 is an asymmetric monomer containing two aromatic and two polar features held in a distinct shape by a peptide-like γ-turn (Figure 2). This minimum pharmacophore was consistent with results reported by Pharmasset, Inc. scientists (now Gilead Sciences, Inc.) regarding the importance of the γ-turn motif in their fluoro-olefin peptidomimetics.[37]

Automated Docking to the Homodimeric Receptor Models

Automated docking of compound 1 was performed as described in Experimental Methods resulting in identification of 10 low-scoring binding poses for each potential receptor model for a total 200 poses. Poses with higher interaction scores appeared to bind peripheral sites along the AH, while the poses with the lowest energy interactions tended to occur closer to the dimeric core. Since the correct binding pose does not always have the absolute lowest docking score,[38] we further ranked the lowest scoring receptor/pose sets from both the fixed and flexible NS5A models for match to the minimum pharmacophore, similar to a procedure recently done for optimized kinase receptor modeling.[39] The top two results were distinct, low-energy, asymmetric binding orientations (mode-I and mode-II) that are both consistent with our library-derived pharmacophore (Figure 3). Each binding mode involves the symmetric caps of compound 1 binding to two distinctly different sites associated with residues 93 and 31 shown in space-filling representation. In mode-1, γ-turn aligned rings A′, B′, and C′ of compound 1 match the pharmacophore and orient the flexible carbamate feature of D′ into a central site at the protein dimer core with potential for H-bond bridging between residues Y93 of either monomer (site 1). The second cap of compound 1 is packed against a complementary steric surface of L31 at the Y93 dimer interface in this receptor conformation. The biphenyl linker lies within a hydrophobic cleft formed above P35 and P32 at the extended PxxPxxP dimer interface. In mode-II, rings A, B, and C of compound 1 changed conformation to match the pharmacophore γ-turn and placed the D carbamate within a site between residues Y93 and L31 of opposite chains that is revealed by concerted hinge-like movements of the PxxPxxP linkers and AH of each chain relative to D-Ia (site 2). Specific interactions of the cap within site 1 change because of the different conformation and orientation of mode-II.

Figure 3

Development of structure-based models for evaluation of activity relations. Best-ranked two binding modes for 1 are at the AH/D-Ia dimer interface. Mode-I: The monomeric pharmacophore features of Figure 2 are inserted into a deep pocket between A-chain Y93 (gold) and B-chain Y93 (blue) at the core of the NS5A-D-I homodimer. The remainder of compound 1 binds against a complementary surface of L31 at the AH interface but is partially exposed and thought to be of lower affinity. Mode-II: The monomeric pharmacophore features fit tightly within a cleft between Y93 and L31 of opposite monomers resulting from a hingelike movement of P35 near the dimer core that shifts the PxxPxxP linker motif.

N-Term Orientation and Asymmetric Binding Offer Shared Role for Positions 93 and 31 in Drug Resistance

Supporting Information Figure S-3 provides a more detailed view of the two sites involved in compound 1 binding. Site 1 is located at the core interface defined by positions of residues Y93 and F36 from the Tellinghuisen et al. dimer structure and the inward directed positioning of N-terminal residues of AH modeled using the SNX5 template. AH folds across residues L31 and L28 as spacers above the mobile PxxPxxP motif. Concerted, asymmetric movements of the AH-PxxPxxP linker region (residues P29–P35) relative to the compliment (residues Y93–P97) region of D-Ia create the second site between 93 and 31 of opposite chains and modify N-term interactions at site 1. A movie of the proposed concerted movement between the two binding modes is provided as Supporting Information. Accordingly, mutations at positions 93 and 31 can be coupled to drug binding changes in either site. Since both modes had similar predicted Autodock binding scores, we evaluated each model for consistency with potential biological mechanisms of action.

Hypothetical Two-Site Model of Asymmetric Drug Binding

Our analysis of binding across the receptor models resulted in a recently published two-site hypothesis that may be useful in understanding the broad chemical space of this class and subcellular redistribution associated with their high activity.[2] Our structural hypothesis is consistent with modes of action (MOA) reported by Targett-Adams et al. that suggest NS5A-DAA perturbs the function of new replicon complexes rather than acting on preformed complexes, and causes redistribution of NS5A to lipid droplets (LD) in cells.[6] The potential membrane alignments using three conformational states extracted from the 20 NS5A-D-I receptor models are illustrated in Figure 4, Figure S-1, III. Mode-I places primary interaction of compound 1 at the core site between Y93 of each monomeric subunit (A-Y93-B) with the second cap against linker residue L31 at the extended dimeric receptor interface formed by AH (Figure 4a). Drug binding in this endoplasmic reticulum (ER) membrane-aligned mode may inhibit cofactor complexes and reversible invagination of NS5A needed to form the replication complex (RC). A conformation of AH that would be most consistent with preformed RC within the ER is shown in Figure 4b. The docking of compound 1 to this “internally curved” RC receptor model gave no binding to either site, forming only low-affinity interactions. Conversely, the mode-II binding model to the hinge site 2 (Figure 4c) closely fits an “externally curved” membrane that could be associated with binding lipid droplets (LD) and vesicles budding into cytosol. Therefore, mode-II binding can provide a 2-fold mechanism to account for the unique picomolar activity of this class: (1) Drug binding in mode-II provides a high-affinity wedge inhibiting the intramolecular protein–protein interactions needed for NS5A to form replication complexes. (2) Drug binding between the mobile PxxPxxP linker, AH, and core stabilizes a protein conformation associated with irreversible sequestering of NS5A in LD related structures. Accordingly, we believe that the unique high potency of this class may be due to a structural mechanism of driving NS5A receptors away from the ER. This two-stage structural MOA may help explain two modes of action for these compounds and their associated kinetics of RC formation and virus particle release observed in both cell-based and clinical studies.[9−11,40]

Figure 4

Proposed structural mechanism of compound 1 activity. Drug (cyan) binds simultaneously to two asymmetric sites at the NS5A D-I A/B dimer membrane interface (blue/gold). (a) The drug binds dimer at ER-membrane surface with primary interaction at A-Y93-B core binding site and a lower affinity site near Leu31 formed by AH in plane with the ER membrane. (b) NS5A binding with RNA, NS5B, and other proteins induces conformational change of AH associated with membrane bending toward replication complex (RC) formation. Replication complex and membranous web formation significantly lowers affinity for drug binding. (c) Conformational movement of AH exposes a drug site at a hinge region between 93 and 31 of the different subunits. Binding to this hinge site locks NS5A D-I into a bent conformation conducive to lipid droplet formation and release to cytosol and is thought to impair assembly of other viral oligomers. Mode-II binding is thought to be associated with the high activity of this class.

Mapping of SAR “Activity Cliffs” to Mode-I and Mode-II

To further evaluate these theoretical models, we tested them using our chemical library of over 67 structurally related analogs of compound 1 published previously.[31−33] We used chemotype clustering and “activity cliff” assessment to reveal trends and summarize results across all the compounds in our three scaffold series.[41] We defined an activity cliff as compound pairs having shared scaffolds with only a single group or atom change that causes >10-fold change in EC50. The maximal cliffs within each of our three scaffold series are summarized in Tables 2–4 and were fit within each binding mode to differentiate high-activity (pM) from low-activity (nM) and inactive (μM) analogs (Figures 5–7).

Table 2

Series I Activity Cliff SAR: Rings A/B Substituted Monomeric Analogs

Table 4

Series III Activity Cliff SAR of C Ring Dimer Analogs

Figure 5

Series I monomeric analogs: mode-I and mode-II binding of activity cliffs. Chains A and B of the NS5A D-I dimer are represented as gold and blue ribbons, respectively. Interaction surface of residues within 4.5 Å of the ligand is colored for hydrophobicity (green) and hydrophilicity (magenta). (A) Binding mode-I of compound 5, with γ-turn at core 93–93 dimer interface, is consistent with placement of monomeric pharmacophore hydrophobic features (gold spheres) of rings A and B in hydrophobic groove formed between F36 and Y93 of different subunits at the dimer interface. A-ring azido substitution (blue atoms) provides additional favorable packing with Y93. Compound 5 was selected for mutation Y93H in our HCV-1b replicons. (B) Binding mode-II places γ-turn at 93–31 dimer interface with ring C packed against proline (P29) of chain A (gold ribbons) while still allowing favorable interactions with Y93 for the A-ring azido. Structural rationale for 10-fold potency increase of compound 6 may be due to conformational stabilization and placement of ring B substitution (Br, red atom) in a halogen-binding pocket at L31/AH interface. Consistent with this interaction model, compound 6 was selected for mutation L31F in our HCV-1b replicons (Table 2).

Figure 7

Series III activity cliff: stereochemical effect of C-ring analogs of compound 1. Chains A and B of the NS5A D-I dimer are represented as gold and blue ribbons, respectively. Interaction surface of residues within 4.5 Å of ligand is colored for hydrophobicity (green) and hydrophilicity (magenta). Left panel: Both compounds 9( and 10( are predicted to bind in mode-II with R4 group, azido, or carboxymethyltriazole, fitting near a complementary surface of chain B, consistent with observed high potency. Right panel: (C) While the model predicts compound 9( can bind in mode-II with C-ring azido favorably fitting against the chain A receptor surface, the bulky carboxymethyltriazole of compound 10( is predicted to clash with the receptor at P29, consistent with complete loss of activity observed for this epimer (Table 4).

Series I: SAR and Binding of A and B Ring Substituted Monomer Analogs

The broad SAR of 36 analogs in our series I that systematically explored nondimeric modifications of compound 1 has been discussed previously, and the largest activity cliffs are summarized in Table 2.[31] The monomeric scaffold compound 4 displayed no activity at micromolar concentrations in our replicon assay. However, a large activity increase was observed for phenyl or azido substitutions selectively placed at the para position of ring A (R1) while ortho- or meta-substitutions with the same groups remained inactive.[31] Our derivative compound 5 is representative of this large cliff showing an activity increase to 8 nM for the para azido from the inactive parent. Compound 6 highlights a smaller cliff associated with halogenation of ring B. While all halogen substitutions showed modest effect, Br substitution at position R2 of compound 6 caused a >10-fold increase in potency compared to compound 5, exhibiting picomolar inhibition of genotype 1b replicons. Resistance selection of Y93H and L31F for compounds 5 and 6, respectively, suggests that these monomeric fragments both act on the same NS5A target site as dimeric compound 1 (Table 2). Mode-I and mode-II binding of compounds 5 and 6 illustrate a structural rational for the observed activity cliffs and resistance selection (Figure 5).

Mode-I binding of compound 5, shown in Figure 5a, places its proline–valine–carbamate cap (rings C and D) within site 1 formed by residues Y93 and F36 of the two different protein subunits, chains A and B, between the Zn2+ binding sites at the dimer interface. The drug’s proline-like ring C stabilizes an internal γ-turn that closely matches our pharmacophore. It also positions compound 5’s C and D amide groups for hydrogen-bond donor/acceptor pairing (blue spheres) with the Y93 hydroxyls of both chains. The compound’s hydrophobic D valine packs against the complementary surface of NS5A chain B residue B-F36 (green), while its ring C γ-turn further aligns hydrophobic rings B and A (gold spheres) within a matched hydrophobic channel at the interface of the N-terminus of chain A and Y93 of the protein dimer’s chain B. Compound 5’s ring B (gold spheres) packs between hydrophobic residues A-F36 and A-P35 (green surface) to place the compound’s ring A between chain A-P35 and chain B-Y93. The para-azido substitution of ring A is aligned for stacking with residue B-Y93 and complementary interactions with AH and linker residues of chain A. Accordingly, the azido groups are predicted to increase affinity for compound 4 only in the para position, consistent with the experimentally observed cliffs.[31]

Mode-II alignment of compound 6 (Figure 5b) places its proline–valine–carbamate cap into the second complementary site within the cleft of the PxxPxxP AH-linker “hinge” region revealed by our flexible modeling. The hydrogen-bond pharmacophore features (blue spheres) are oriented by the γ-turn to form strong buried hydrogen bonds between chain B-T95 with donor/acceptors of the drug’s rings B and C, while the carbamate moiety (portion D) interacts with the solvent exposed arginine, A-R30. The proline-like ring C of compound 6 packs against the proline P29 receptor surface, aligning its hydrophobic rings B and A between the receptor’s chain A linker residues A-P32, A-P35, and chain B-Y93 at the dimer interface. This orientation also places the drug’s para-azido substitution of ring A for stacking with residue B-Y93. Additionally, it reveals a complementary pocket for compound 6’s bromine (ring B, red atom) at the L31 surface that is predicted to increase relative affinity for this mode. The increased affinity for this hinge-binding mode may drive the 10-fold increased potency of this substitution in line with our asymmetric two-site mechanism and is consistent with rapid selection of the L31F mutation for compound 6 compared to Y93H for compound 5 (Table 2).

Series II: SAR and Binding of A and B Ring Substituted Dimer Analogs

Series II explored 12 analogs with scaffold modifications that included replacement of the two phenyl rings with shorter, nonaromatic linkers that resulted in >1000-fold drop in potency, while extensions of the aromatic linker A–A′ by up to five rings retained picomolar potency (not shown).[33] These data further support the asymmetric binding mode-II model associated with high potency and the mechanism suggested by Figure 4c. The ∼1000-fold activity cliff presented in Table 3 for compounds 7 and 8 with a bent R3 linker substitution was also consistent with the mode-II binding model as illustrated in Figure 6. The modeling suggested that halogen substitution at B-ring of compound 8 increased the potency of compound 7 through improved binding affinity for mode-II in two ways. First the binding mode placed the ring B bromine of compound 8 within the halogen pocket near L31, which is similar to monomeric compound 6, and the ring B′ bromine substitution fit another halogen complementing pocket at the Y93 core interface. Second, conformational searching of both analogs 7 and 8 revealed that bromine substitutions of the B rings restricted the floppy compound 7 to a deep energy minimum near the conformation shown for compound 8 bound in mode-II (Figure 6). Therefore, our model suggests the large potency increase for compound 8 may be due to favorable changes in both enthalpy and entropy associated with binding in mode-II.

Table 3

Series II Activity Cliff SAR of A′–A/B Linker Modified Dimers

Figure 6

Series II A′-A/B linker modified dimer activity cliff: mode-II hinge-binding. Chains A and B of the NS5A D-I dimer are represented as gold and blue ribbons, respectively. Interaction surface of residues within 4.5 Å of ligand is colored for hydrophobicity (green) and hydrophilicity (magenta). Both compounds 7 and 8 can only partially bind in mode-I. However, mode-II fitting of compound 8 (shown) provides halogen binding pockets for both R2 and R2′ Br substitutions (red atoms of rings B) via the “bent” R3 linker and stabilizes this minimum energy conformation consistent with its higher activity (Table 3).

Series III: SAR and Binding of C Ring Substituted Dimer Analogs

For series III including 19 analogs with modifications of rings B and C, the largest activity cliff occurred on ring C.[32] Stereoselective azido-group substitution of the proline-like ring C did not significantly change the inhibitory activity of compounds 9( or 9( from that of compound 1. However, the larger 1,4 triazole groups of compounds 10( or 10( were only tolerated as (R) conformation, suggesting stereochemical restraints of the binding pocket (Table 4).

Our mode-II binding models of these compounds showed that substitution of ring C with an azido-group can be comfortably accommodated for compound 9() and a potential hydrogen-bond to the backbone NH of P29 may occur for the closely packing compound 9() (Figure 7, top panels). Consistent with these results, there is room for the bulky triazole substitution of compound 10() to fit and maintain high potency; however, compound 10() clashes with the surface of P29 leading to a complete loss of activity for that isomer.

Modeling Other Clinical Candidates

To validate the relevance of our models using compounds outside our training set, we applied the identical force-field based fitting method used for our activity cliff analysis to clinical compounds 2 (GSK-2336805) and 3 (GS-5885).

GSK’s compound 2 fits the hinge-binding mode-II model like 1, with the spiroketal replacement of ring C receiving an additional H-bond from the backbone NH of R30 similar to our compound 9() and complementary packing with P29 (Figure 8). While our mechanism does not predict significant activity differences between these compounds in WT system, the spiroketal H-bond interaction may be the source of the improved L31V and Y93H activities recently reported for compound 2 (3-fold resistance against each mutation for 2 compared to 25- and 64-fold, respectively, for 1; activities against the double mutant were not reported).[15] A potential mechanistic role for L31 positioning the membrane bound AH relative to the core Y93 site is illustrated in Figure 8. Conserved Trp (W) residues (4, 9, 11, 22) of AH align within the cytosolic membrane surface and L31 packs against residues W11 and D10 as a spacer between AH and the PxxPxxP linker. Change of L31 to V or F, or Y93 to H, could slightly change the orientation of AH at both the 31/93 and 93/93 dimer interfaces. This provides a structural relation for the combined effect of mutations at 31 and 93 and suggests the strong H-bond from the backbone to the spiroketal of compound 2 may allow the drug to remain bound to the mutants with greater affinity than compound 1.

Figure 8

Mode-II: hinge binding model of clinical compound 2 (GSK-2336805). Translucent van der Waal surface of residues 31 and 93 are colored for hydrophobicity (green) and hydrophilicity (magenta). GSK’s compound 2 binds at the cytosolic membrane interface of amphiphillic helix (AH). Spiroketal replacement of ring C affords an additional strong H-bond with the backbone NH of R30 similar to our 9( and close packing with P29. These enhanced interactions may increase the barrier to 31 or 93 resistance relative to compound 1 by securing position between chains A and B of D-I dimer represented as gold and blue ribbons, respectively. Distance labels between residues P29 and P97 at the asymmetric mobile AH-D-I linker interfaces are shown in green.

Gilead Sciences compound 3’s cyclopropyl addition to ring C proline can also be accommodated similarly to GSK’s compound 2 but without the backbone H-bond (Figure 9). The remainder of the tricyclic linked, pseudodimeric compound 3 finds new interactions within our asymmetric model. For example, mode-II binding of the tricyclic benzimidazole–difluorofluorene–imidazole core places the fluorines into a complementing pocket formed by the methylene of aspartic acid residue D10 on AH. Consistent with recently published SAR for the compound 3, the dihalogen substitution appears optimal for the tricyclic bridge of the A–A′ rings, but larger groups such as O or Me can also be accommodated.[14] The activity of the asymmetric B′ and C′ group substitutions can also be explained with B′ acting as H-bond acceptor for residue R6 of AH and orienting the C′ and D′ groups similarly to compound 1 within the dimer core between Y93 of each chain. This interaction positions the bulky bicyclic C′ group of Gilead’s compound 3 in a solvent exposed space between the N-term residues, S3, flanking the binding site. This alignment of the drug relative to residues 28, 30, 31, and 93 provides structural context to model resistance to compound 3 observed in recent clinical trials.[42]

Figure 9

Mode-II: hinge binding model of clinical compound 3 (GS-5885). Translucent van der Waal surfaces of residues 3, 28, 31, and 93 are colored for hydrophobicity (green) and hydrophilicity (magenta). Chains A and B of D-I dimer are represented as gold and blue ribbons, respectively. Gilead Sciences compound 3’s cyclopropyl ring C is accommodated like compounds 1 and 2; the tricyclic A–A′ analog places halogen within a pocket of the AH. Pseudodimer B′–C′–D′ fits within this asymmetric model, consistent with its high activity and recently reported SAR.

Comparison to Existing Models

A recent study of potent NS5A-DAA binding suggests a major mode-of-action for these drugs is blocking formation of replication complexes within the membrane surface. Their asymmetric models of binding to the “clamlike” 1ZH1 dimer form near residues 93 and 54 are very close to ours, and the authors suggest binding to this “membrane-proximal surface” may disturb “positioning and/or folding of the N-terminal linker segment connecting DI with AH” that is consistent with our proposed model.[11] A study from BMS also finds the “best fit” of their direct photolabeling data using the 1ZH1 dimer form. While their symmetric binding model is similar to some of ours shown in Figure 3, docking the most potent stilbene based analog used in their experiments with our flexible receptor method finds a best fit to the same “hinge binding” receptor as compound 1 placing the azido moiety near residues 30 and 29, consistent with their labeling results.[20]Supporting Information Figure S-4 illustrates binding of the stilbene scaffold and three other recently reported active chemotypes docked and ranked using our method. Although researchers at GSK developed a different model based upon the alternative 3FQQ dimer form, our result for their lead compound illustrated in Figure 8 is consistent with the SAR and resistance data recently reported.[15]

Conclusions

We have developed a novel semiautomated method for production of in silico models to predict conformations of the full length NS5A D-I protein receptor suitable for automated docking and SAR evaluation. Our method incorporates both a novel conformation of the conserved PxxPxxP AH-D-I linker motif and a core-directed N-terminal AH fold based upon a structurally homologous membrane binding protein (SNX5) as empirical constraints while modeling flexible biological space. We developed a minimum pharmacophore to probe and rank docking models associated with drug activities. These results suggested two novel asymmetric binding modes that are consistent with much of the recent scientific literature for NS5A-DAA mechanism(s) of action including direct photolabeling by NS5A-DAA,[20] inhibition of RC formation at the ER, and sequestering NS5A in lipid droplets to inhibit virus particle formation/release by drug binding at the membrane-bound AH-D-I interface.[11] We used “activity cliffs” from our analog SAR library and other clinical compounds as validation of the binding modes. These models provide explanation for the high activity of this class and the specific effects of selected drug-resistant mutations. These methods and models are presented as tools toward understanding the unusual activity of this class with potential to aid development of NS5A-DAA with increased barrier to resistance and improve treatment of HCV infection. (3D coordinates for NS5A-Gt1b mode-I and mode-II models used for this study are available upon request).

Experimental Methods

Analysis of Experimental Structures

All experimental coordinates for NS5A D-I subdomains, used as references for full-length model building, were downloaded from the RCSB Protein Data Bank (PDB) at www.pdb.org.[43] NMR ensembles and minimized average structures of genotype 1a N-terminal residues 1–31 that encode the membrane binding the amphipathic helix (AH) (residues 5–25) were downloaded as 1R7C, 1R7D, 1R7E, 1R7F, 1R7G. X-ray crystal structures of subdomains D-Ia/D-Ib from genotype 1b were downloaded in two dimeric forms, 1ZH1 (residues 36–198) and 3FQQ (residues 31–191). Structures of 1ZH1 chain A and 3FQQ chain B were aligned using the Matchmaker function in Chimera[44] with these parameters. Chain pairing: bb, Needleman–Wunsch using BLOSUM-62. ss fraction: 0.3, gap open (HH/SS/other) 18/18/6, extend 1. ss matrix: (O, S) −6, (H, O) −6, (H, H) 6, (S, S) 6, (H, S) −9, (O, O) 4. Iteration cutoff: 2. The rmsd between the 154 atom pairs matched was 0.538 Å. Of note, the most divergent residues observed from the aligned monomers were between F36 and F37 with α-carbon distances of 6.2 and 9.2 Å, respectively.

Structure Based Sequence Alignments

Sequences for HCV genotype 1a and 1b NS5a were downloaded from Viral Bioinformatics Resource Center (VBRC) (http://www.hcvdb.org/). Sequences of each unique structural subdomain were aligned for comparison with full-length sequences of NS5A 1a and 1b in UCF Chimera.[44] The sequence of file 1R7G is representative all NMR-based AH-subdomain structures evaluated. Domain motifs, D-I (residues 1–213), D-II (residues 250–342), and D-III (residues 356–447), as first defined by Tellinghuisen et al.,[26] were labeled, and positions of drug resistant mutations 31 and 93 were mapped to sequences (Figure 1). We defined residues 28–36 (containing the L31 residue) to be the flexible AH-D-Ia linker (Figure 1, orange box).

Flexible AH-D-Ia Linker Modeling Workflow

D-Ia/D-Ib dimers, 1ZH1 and 3FQQ, were aligned with the Matchmaker function in Chimera.[44] AH models were manually aligned to each dimeric structure in Maestro, version 8.1 (Schrödinger, LLC, New York, NY); however, no biochemically reasonable solution was found to link experimental AH models consistent with NMR and drug resistance data. De novo automated building of residues 1–36 from D-Ia/D-Ib template was also unsuccessful. Subsequently, we narrowed our conformational sampling to the AH-D-Ia linker containing a conserved PxxPxxP motif (Figure S-1, Ia).

Modeller, version 9.10,[34] was run with default settings to extend the octapeptide from each experimental D-Ia/D-Ib dimer template (1ZH1 and 3FQQ) (Figure S-1, Ib). The 20 lowest energy, unique conformers were retained for each dimer model. Ramachandran plots were generated in Maestro and Chimera[44] to check protein health. Intradimer and interdimer distances between 31/93 α carbons were calculated and ranked using a custom script for Chimera (Figure S-1, Ic). The conformation with the lowest distance (9.8 Å) was located on one 1ZH1 monomer and used to generate a symmetrically extended dimer for next stage AH fitting (Figure S-1, Id).

AH-Linker Homology Based on PxxPxxP

We performed a search for experimentally solved proteins having PxxPxxP motifs to evaluate potential relevance of our novel linker conformation. The phox (PX) domain of sorting nexin 5 (SNX5) was identified, and associated coordinates of structures 3HPB and 3HPC were downloaded from PDB.[35] Each set contains the motif at P91, P94, and P97. Only small differences exist between the two structures, solved under different conditions, and we chose to use both of them. Each structure was fit to one PxxPxxP motif, P29, P32, P35, of the symmetrical, linker extended homodimer generated in previous workflow (Figure S-1, IIa). Three atoms of each proline residue Ca, Cg, N, were fit via command line “match” function in UCSF Chimera, version 1.6.2. 3HPB (green) residues P91, P94, P97 were fit to P29, P32, P35 of chain B and of our homology model with rmsd of 1.5 Å between nine atom pairs. 3HPC (salmon) residues P91, P94, P97 were fit to P29, P32, P35 of chain A of our homology model with rmsd of 1.5 Å between nine atom pairs. Both fit 3HPB and 3HBC structures were trimmed back to residues His69-Phe99 as shown (Figure S-1, IIb).

NMR Based AH Fitting on SNX5 Template

1R7E.pdb, a minimized average NMR structure of the membrane anchor domain (1–31) (AH) of the NS5A was fit to each SNX5 truncate using Matchmaker in Chimera [chain pairing: bb; Needleman–Wunsch using BLOSUM-62. ss fraction: 0.3; gap open (HH/SS/other) 18/18/6, extend 1. ss matrix: (O, S) −6, (H, O) −6, (H, H) 6, (S, S) 6, (H, S) −9, (O, O) 4. Iteration cutoff: 2]. Notably, 1R7E match to 3HPC rmsd between 14 atom pairs is 0.939 Å, and 1R7E match to 3HPB rmsd between 13 atom pairs is 0.686 Å (Figure S-1, IIc).

Membrane Alignment of AH

Alignment of the fit AH was evaluated in the context of a POPC model popc128a.pdb downloaded from the Tieleman laboratory http://people.ucalgary.ca/~tieleman/download.html (Figure S-1, IId).[45]

Sequence Based Building of Full-Length NS5a D-I

The aligned, trimmed, SNX5 dimer was used as structural template for homology modeling in UCSF Modeller, version 9.10.[34] The sequence of HCV genotype 1b NS5A, amino acids 1–198, from VBRC was used to build a full-length model of NS5a D-I dimer. The single contiguous result was then used as the structural template to explore potential flexibility under two classes of conditions: (1) Coordinates of 1–198 from the previous full-length Gt1b model were used as template to model potential “fixed” conformations. (2) Coordinates 1–38 of the Gt1b model and 36–198 of 1ZH1 were both input as templates allowing AH to move around residues 36 and 37 as a “flexible” motif. Ten dimer conformers were generated for each class suitable for automated docking (Figure S-1, III).

Chemical Library Development

Sixty-seven analogs were synthesized in three congeneric series and tested as described previously (Figure S-2).[31−33] Briefly, the three series were tested for activity in HCV genotype 1b Huh-7 subgenomic replicons with parallel toxicity in Huh-7 and other cell panels including primary lymphocytes. Chemical structures and activity data were joined to create a compound library for SAR evaluation using a custom workflow in Pipeline Pilot (Accelrys Corp, San Diego, CA). Joined data were imported into MOE 2011.10 (Chemical Computing Group, Montreal, Quebec, Canada) and analyzed with SAReport module.

The full set was filtered for activity of <0.99 nM, resulting in 25 compounds with at least one representative from each structural set shown in Figure 2 and discussed in results section. Conformations for the compounds were generated in MOE 2011.10 (Chemical Computing Group, Montreal, Quebec, Canada) using LowModeMD with default MMFF94x force field and rejection and iteration limits of 100 and 50, respectively. The rms gradient cutoff of 0.5, MM iteration limit of 10, rmsd limit of 1, energy window 7, and conformation limit of 30 were used. MOE pharmacophore elucidation was performed with the clustered conformations using a unified scheme and default features.

Automated Docking and Ranking of Daclatasvir Binding

A large grid size of 36 Å × 76 Å × 36 Å was defined in Autodock tools (http://autodock.scripps.edu/resources/adt) to explore potential for drug interaction with the entire AH-D-Ia surface of all 20 conformations illustrated in Figure 3 and discussed in results section. Autodock Vina (http://vina.scripps.edu) was used to place and score daclatasvir 1 within the grid, and the 10 lowest energy-binding modes were saved for each receptor conformation.[46] Receptors were ranked by interaction energy and the lowest scoring receptors were run again using a smaller grid of 26 Å × 76 Å × 20 Å to optimize lowest energy complexes. The binding poses of compound 1 with lowest-energy Autodock scores were further ranked by manual fitting to the best pharmacophore derived above and used for SAR evaluation.

SAR Evaluation of Activity Cliffs

The 67 compounds from our chemical library were clustered by synthetic scaffold, and each cluster was evaluated for “activity cliffs” similar to our earlier work.[47,48] We defined activity cliffs as compound pairs having shared scaffolds with only a single group or atom change that causes >10-fold change in EC50. For each of the three analog classes, those pairs with most similar structures and highest potency differences are presented in Tables 2–4 and were modeled within the receptor for SAR analysis. Top ranking Autodock model complexes were read as PDB into MOE 2013.08 (Chemical Computing Group, Montreal, Quebec, Canada). Systems were checked for protein health and atom typing before addition of hydrogens using the Protonate 3D function. Ligand and protein were subjected to stepwise tethered minimization using the MMFF94x force field. The Rotamer Explorer function was used to test alternative side chain rotomers and Protonate 3D was reapplied before each new test of minimization as described previously.[49] Novel analogs were built from the monomeric binding of mode-I (93–93 site 1) and mode-II (93–31 site 2) and tested for site shape complementarity and potential electrostatic interactions using a similar stepwise minimization scheme. Final minimizations for each were performed for all atoms of ligand and residues within 4.5 Å using no tethering and gradient of 0.1 rms kcal mol–1 Å–2. The activity cliffs used for presentation were found representative of all SAR observed in the larger sets.

Resistance Selection

Compounds 1, 5, and 6 were submitted for resistance selection in HCV Huh-7 subgenomic replicons containing 10 000 cells/well in six-well plates with 1 mg/mL G418 in selection media. Titrations were performed with selecting agents at 0.5, 1, 2, 5, and 10 times the EC90 concentration. RNA was isolated approximately 2 months after selection from single clones or remaining population and sequenced. Results are presented in Table 2.