Low susceptibility of asunaprevir towards R155K and D168A point mutations in HCV NS3/4A protease: A molecular dynamics simulation

Abstract

Hepatitis C has become an important health problem that requires expensive treatment and leads to liver tumorigenesis. Hepatitis C virus (HCV), which is the main cause of hepatitis C, has a high mutation rate due to the lack of proofreading activity of the RNA polymerase enzyme. The NS3/4A serine protease is an important target for anti-HCV drug discovery and development because of its crucial role in the cleavage of the polypeptides involved in viral replication. In the present study, all-atom molecular dynamics simulation was performed to elucidate the effect of the single point mutations R155K and D168A in the HCV genotype 1 NS3/4A protease on the structural dynamics, molecular interactions and susceptibility of asunaprevir (ASV), a second-generation NS3/4A protease inhibitor. Principal component analysis indi- cated that these two mutations converted the direction of motion of residues 123, 155 and 168 in the binding pocket to significantly point outwards from ASV, resulting in a loss of the hydrogen bond network of residues R123$$$R155$$$D168. The free energy calculations based on different semiempirical QM/MM-GBSA methods revealed that the binding affinity of ASV with the two mutant forms of the NS3/ 4A protease was significantly decreased in the order of wild-type < R155K < D168A. This work provided useful structural information regarding the atomistic understanding of acquired drug resistance against ASV caused by the R155K and D168A mutations.

1. Introduction

Hepatitis C virus (HCV) is a small (50e80 nm in size) enveloped RNA virus belonging to the Hepacivirus genus of the Flaviviridae family [1]. The HCV particle consists of a nucleocapsid, which contains the positive single-stranded RNA genome covered by a host cell-derived lipid envelope [2]. Currently, HCV infects over 80.2 million people and long-term HCV-infected patients are under an increased risk of developing liver diseases, such as cirrhosis and hepatocellular carcinoma. Genotype 1 is the most common HCV genotype worldwide, accounting for 46% of all HCV infections fol- lowed by 30% from genotype 3. The highest diversity of genotypes is observed in China and South-East Asia, while only genotype 4 is found in Egypt and Mongolia [3e5]. The currently used anti-HCV agents are targeted inhibitors against the important NS3/NS4A protease, NS5B polymerase and NS5A viral proteins [6]. For the treatment of HCV patients, the combination of PEGylated-inter- feron-a and ribavirin is extremely expensive but only effective for
~50% of HCV patients [7,8]. Moreover, this therapy causes several adverse side effects, including alopecia, rash/itching, nausea, thrombocytopenia and anaemia.

In the HCV life cycle, after protein translation, the large poly- peptide encoding the nonstructural proteins is processed by the viral NS2/3 and NS3/4A proteases. The NS3/4A protease cleaves the NS3A-NS4A, NS4A-NS4B, NS4B-NS5A and NS5A-NS4B junctions and so is essential in the viral replication process. Consequentially, the NS3/4A protease has become one of the attractive targets for anti-HCV drug design and development [9]. Importantly, the NS3 proteolytic activity and the acceleration of the cleavage rate are associated with the NS4A cofactor [10]. The HCV therapy has been improved by the use of boceprevir [11e13] and telaprevir [14,15], which specifically inhibit the NS3/4A protease enzyme. Although this treatment leads to an increased rate of sustained virologic response, such therapeutic drugs have a low success rate [16]. Recent drugs, such as simeprevir [17,18] and vaniprevir [19,20], were approved for the treatment of HCV-infected patients in 2013 and 2014, respectively. In addition, other inhibitors are currently used in clinical studies such as danoprevir [21,22], grazoprevir [23], faldaprevir [24], vedroprevir [8] and asunaprevir (ASV; [25,26]). Because HCV lacks a proofreading activity of its RNA polymerase enzyme, high mutation rates in response to the treatments have been widely reported, resulting in resistance to drugs/inhibitors [27e30].

Asunaprevir (ASV) is a second-generation NS3/4A protease inhibitor that was approved for use in combination with daclatasvir in Japan for the treatment of HCV genotype 1-infected patients [31]. However, ASV has not yet been approved in the USA and Europe, where it is currently in phase III clinical trials [32]. This compound shows a strong antiviral activity against genotypes 1 and 4 of HCV [25,26]. From in vitro studies, the acquired resistance against ASV caused by the two single point mutations of R155K and D168A in the NS3/4A protease results in a lower susceptibility to ASV, with an
~21- and 23-fold reduced half maximal effective concentration (EC50), respectively [28]. These two mutated residues are located near the protease active site (Fig. 1) and were also found to reduce the binding affinities of some other anti-HCV inhibitors, such as danoprevir [21,22], simeprevir [17,18] and grazoprevir [23], by changing the conformation of the active site [33e35]. For instance, the absence of a Nε, which is the nitrogen atom at epsilon position, of residue K155 prevented the interaction between residues 155 and 168, which disrupted the interaction of ASV with the target protein [36].

To investigate the effect of such mutations on the ASV binding to NS3/4A, in terms of the inhibitor-protein interactions and binding efficiency towards the HCV NS3/4A protease, all-atom molecular dynamics simulations (MDSs) and binding free energy calculations based on QM/MM-GBSA and MM/GB(PB)SA methods were applied on the ASV in complex with the wild type (WT) and the R155K and D168A mutations of the HCV genotype 1 NS3/4A protease (Fig. 1). In addition, the crucial motion of amino acids located in the active site of the apo (APO) and complex (CPX) forms was characterized using principal component analysis (PCA).

2. Methods

2.1. System preparation and ligand optimization

The three-dimensional structure of ASV in complex with the WT (PDB: 4WF8 and R155K mutant (PDB: 4WH6)) NS3/4A forms [36] was obtained from the protein data bank. The enzyme variant containing the D168A mutation was built by changing residue 168 from aspartate (D) to alanine (A) using the Discovery Studio 2.5 Accelrys Inc. [37]. The apo proteins were generated in all variants by removing the ligand from the complex structure. The protonation state of all ionizable residues was characterized using PROPKA3.1 [38] at pH 7.0. Note that, the cysteine (C) residue that is coordinated with Zn2+ was set as deprotonated state (CYM). The missing hydrogen atoms were added using the LeaP module [39] in AMBER14 [40]. The electrostatic potential charges of the ligand were computed at the HF/6-31(d) level of theory, while the restrained ESP charges and corresponding parameters of the ligands were generated using the antechamber and parmchk modules, respectively, in AMBER14 [41e43]. The AMBER ff03.r1 and GAFF force fields were applied for protein and ligand, respec- tively. Chloride ions were added for neutralizing the simulated systems. Each complex was solvated in a 80.9 × 69.5 × 77.9 Å3 rectangular box of TIP3P water [44]. The distance between protein and solvation box edge was 10 Å. Prior to perform the MDSs, the added hydrogen atoms and water molecules were minimized with 1500 steps of steepest descents and conjugated gradient. Lastly, the whole system was minimized using the same procedures.

2.2. Molecular dynamics simulations (MDS)

The MDSs of all systems were performed under a periodic boundary condition in accord with the standard procedures [40]. All covalent bonds involving hydrogen atoms were fixed using the SHAKE algorithm [45]. The cutoff of 10 Å was employed for non- bonded interactions, while the particle mesh Ewald summation method [46] was applied for calculating the long-range electro- static interactions. A Langevin algorithm [47] was used to control the temperature with a collision frequency of 5.0 ps—1 for the first 3.5 ns and subsequently changed to 2.0 ps—1. Each system was heated up to 298 K for 200 ps. Afterwards, all simulations, including the enzyme in APO and in CPX, were treated with the NPT ensemble at this temperature and a pressure of 1 atm using the PMEMD module. All systems were simulated for 100 ns at a time step of 0.002 ps [42]. In each system, 500 snapshots were collected every 1 ns of MDS.

The binding free energy of all systems was performed using the semi-empirical QM/MM-GBSA method compared with MM/GB(PB) SA [48,49]. The inhibitor and the two focused residues (155 and 168) were quantum mechanically treated, while the rest of system was performed using molecular mechanics. The per-residue decomposition free energy (DGresidue) based on the MM/PBSA method was performed to investigate the mutated amino acids (residues 155 and 168) associated with inhibitor binding. Addi- tionally, the important motion of the protein in the WT and mutant (R155K and D168A) systems was characterized by PCA using CPPTRAJ module implemented in AMBER14 [40].

3. Results and discussion

3.1. Stability of the global structure of each simulated model

The all-atom root-mean-square displacement (RMSD) was calculated in order to estimate the stability of each simulated model (WT and the R155K and D168A single point mutations in APO and in CPX with ASV) along the simulation time. As shown in Fig. 2, the RMSD value of each system rapidly increased during the first 10 ns and then fluctuated at ~2.2e2.7 Å in all the studied models until the end of the simulations. Altogether, the MDS sys- tems tended to be stable after 60 ns. In this study, the equilibrated MDS trajectories extracted from the last 15 ns were used for further analysis of the (i) electrostatic network at the mutation site, (ii) protein conformation, (iii) inhibitor-protein hydrogen bond (H- bond) interaction and (iv) effect of mutations on the inhibitor binding affinity.

3.2. Electrostatic network at the mutation site

It has been reported that the electrostatic network formed at residues R123$$$R155$$$D168 plays an important role in the in- hibitor binding [50]. In order to calculate this interaction network throughout the production phase of the simulation, we used the geometric criteria of: (i) a distance between the H-bond donor and acceptor of ≤3.5 Å and (ii) a H-bond angle of ≥120◦. As shown in Fig. 3, a high percentage of H-bond occupation (%HBoc) was detected in the WT system in both the APO and CPX. As expected, the R155K and D168A point mutations importantly decreased the electrostatic network formation compared to that of the WT NS3/ 4A, resulting in a lower binding affinity of ASV towards R155K and D168A NS3/4A. The point mutations decreased not only the elec- trostatic attractions but also the van der Waal (vdW) interactions, indicating that vdW forces also played a role in the drug binding (Fig. 4). The obtained results strongly agreed with previous reports [38,44] as well as our QM/MM-GBSA free energy calculation (as discussed later and Fig. 8).

Fig. 1. (A) Chemical structure of ASV, in which the P10, P2, P3 and P4 substituents are shaded by black, pink, purple, blue and orange, respectively. (B) Three-dimensional structure of the ASV-NS3/4A complex, where the acquired resistance caused by the R155K and D168A point mutations towards the NS3/4A protease is also shown.

3.3. Protein conformation

To understand the effect of the R155K and D168A mutations on the protein conformational changes, PCA was performed [51]. These motions correspond to the correlated vibration or collective motion of groups of atoms. The scatter plots between the first (PC1) and second (PC2) principal components indicated the different distribution of protein conformations.

Fig. 2. RMSD plots of the WT (black), R155K (grey) and D168A (light grey) HCV NS3/4A protease without (APO) and with (CPX) ASV binding for 100 ns of MDS.

The PCA scatter plots (Fig. 5A) revealed that the conformational distributions of the holo (CPX) form of both mutations were dramatically different from those of the WT, implying that the R155K and D168A mutations in NS3/4A promoted protein confor- mational changes in a contrasting manner, which can significantly affect the protein-ligand binding affinity (Fig. 8, as discussed later). The first 15 PC modes of all studied models in the CPX and APO forms showed a percentage of accumulated variance of >50%, indicating that these modes represent the important motions of protein (Fig. 5B).
According to the porcupine plot of the WT system (Fig. 6, top), the direction of motion of residues R123, R155 and D168 moved towards the ligand binding, in a different manner from that in the APO form, suggesting that ligand binding induces the closed form of the protein, in which these three residues are important for ASV recognition. Additionally, the ligand interacted firmly in the pocket site due to the long side chains of two arginine residues. As ex- pected, both R155K and D168A mutations caused the key binding residues 123, 155 and 168 to significantly move outwards from ASV, resulting in a loss of electrostatic attractions and H-bond in- teractions. In summary, these structural transformations could be the reason why ASV showed a low susceptibility towards the R155K and D168A forms of the NS3/4A protease.

3.4. Inhibitor-protein H-bond interaction

One of the important factors involved in the stabilization of protein-ligand complexes is H-bonds. To measure such interactions throughout the overall simulation and equilibration phase (the last 15 ns), H-bond formation was calculated using the two geometric criteria described above, and the results are shown in Fig. 7. Among the five sites of ASV (P10 and P1eP4, Fig. 1A), there were no H-bond interactions between the NS3/4A protease and ASV at the P2 and P4 moieties, in accord with the previous report on faldaprevir and danoprevir recognition [52]. For the WT NS3/4A (black), ASV showed 4 strong H-bonds (%HBoc > 70), one each with G137 and R155 at the P1′ and P1 site, respectively, and two with A157 at the P3 site. With two firm H-bonds detected at the P3 site, the A157 residue played an important role as recognition site for anti-HCV protease inhibitors (Fig. 7A) [42]. Notably, the number of H-bonds formed in the mutated NS3/4A, especially R155K, was lower than those in the WT, suggesting that these mutations attenuate the H- bonds, which is consistent with the loss of electrostatic network formation detected in the R155K and D168A mutations (Fig. 7B).

Fig. 4. Averaged electrostatic and vdW energy contribution from each residue of NS3 protease domain in WT and mutant forms.

Fig. 3. H-bond network of the three important NS3 residues (R123, R155 and D168) in the WT, R155K and D168A systems in both APO and CPX forms.

Fig. 5. (A) The two-dimensional scatter plot of MDS trajectories between the first (PC1) and second (PC2) principal components. (B) The PCA scree plot of quantitative characters of the first 15 modes.

3.5. Effect of the R155K and D168A mutations in NS3/4A on the ASV binding affinity

To evaluate and compare the binding affinity of ASV against the WT and the R155K and D168A mutant strains of the NS3/4A pro- tease, the QM/MM-GBSA binding free energy calculations were performed on 100 snapshots extracted from the last 15 ns of each MDS. Based on this approach, the ASV inhibitor and the NS3/4A residues 155 and 168 were treated by the semiempirical methods (AM1, RM1, PM3 and PM6) compared with MM/GB(PB)SA, while the remainder was described at the MM level. The results of the averaged binding free energies in each system are compared in Fig. 8, which demonstrated that the binding free energy (DGbind) values obtained from all four levels of QM theory are in the same order of WT > R155K > D168A. The R155K and D168A single point mutations resulted in a significant reduction in the binding free energy by ~6 and ~9 kcal/mol, respectively, compared to the WT complex, in reasonable agreement with the experimentally observed 21- and 23-fold changes in the EC50 values, respectively, (McPhee, 2012a). Unlike QM/MM-GBSA method, the free energy values between WT and both mutations obtained from MM/(GB) PBSA calculations were not significantly different. In addition, it should be noted that the calculated binding free energy does not give us an absolutely accurate binding free energy value in com- parison with experimental data, but it evidently shows a relative binding affinity between WT and mutant strains toward ASV. Moreover, EC50 values which refer to concentration of a drug that causes the half maximal effect of an observed effect cannot directly compare to the binding efficiency of drug toward targeted protein. Moreover, the QM/MM-GBSA binding free energy calculations results were similar to previous computational studies, where the binding affinities of narlaprevir and vaniprevir towards the NS3/4A R155K mutant were significantly lower than that of the WT strain [53,54], while the D168 A/V mutants of NS3/4A considerably decreased the susceptibilities of faldaprevir, danoprevir and vani- previr [52,54]. In summary, our QM/MM-GBSA binding free energy calculations not only supported the role of the R155K and D168A mutations towards the weaker binding affinity of ASV, but also strongly predicted the susceptibility of this inhibitor in the same manner as the experimental data.

Fig. 6. Porcupine plot of the PC1 of the studied systems. The arrow head and its length indicate the direction and amplitude of motions, respectively.

Fig. 7. A) Percentage of H-bond occupation over the last 15 ns, and B) the number of H-bond (versus simulation time) of NS3 residues formed with ASV in the WT and the R155K and D168A mutant forms of NS3/4A.

Fig. 8. Binding free energy calculated by the MM/GB(PB)SA and the QM/MM-GBSA method, at the AM1 (black), RM1 (red), PM3 (blue) and PM6 (green) levels of QM theory for the WT and R155K/D168V mutant systems. (For interpretation of the ref- erences to colour in this figure legend, the reader is referred to the Web version of this article.)

4. Conclusions

In this work, 100-ns MDSs were applied on the WT and on the R155K and D168A single point mutations of the NS3/4A protease in the apo form and in complex with ASV, a current drug in phase III clinical trials. According to the PCA, these two mutations dramati- cally converted the direction of motion of the key binding residues (123, 155 and 168) to move outwards from ASV, resulting in a loss of electrostatic attraction and H-bond formation. The QM/MM-GBSA binding free energy calculations at different levels of QM theory (AM1, RM1, PM3 and PM6) suggested that both mutations caused a significant reduction in the binding affinity, ranked in the order of WT < R155K < D168A, which agreed with the observed experi- mental EC50 values. In addition, it can be postulated that other related mutations might have a similar pattern in drug binding and the loss of interactions like in our cases present here. Consequently, our findings can allow us to predict the effect of mutations on the binding affinity at least without simulations. Taken together, the theoretically obtained information can be used as a rational guide for antiviral drug design and development towards the HCV ge- notype 1, which is the most prevalent genotype worldwide.