H274Y's Effect on Oseltamivir Resistance: What Happens Before the

Subscriber access provided by CMU Libraries - http://library.cmich.edu

Article

H274Y’s Effect on Oseltamivir Resistance: What Happens Before the Drug Enters the Binding Site Muhammad Yusuf, Nornisah Mohamed, Suriyati Mohamad, Dusanka Janezic, K.V. Damodaran, and Habibah Wahab J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.5b00331 • Publication Date (Web): 24 Dec 2015 Downloaded from http://pubs.acs.org on December 25, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

H274Y’s Effect on Oseltamivir Resistance: What Happens Before the Drug Enters the Binding Site Muhammad Yusuf,† Nornisah Mohamed†, Suriyati Mohamad,†,‡ Dusanka Janezic,║K.V. Damodaran† and Habibah A. Wahab*,†,┴ †

Pharmaceutical Design and Simulation (PhDS) Laboratory, School of Pharmaceutical Sciences,

Universiti Sains Malaysia, 11800 Minden, Pulau Pinang, Malaysia ‡

School of Biological Sciences, Universiti Sains Malaysia, 11800 Minden, Pulau Pinang,

Malaysia ║

Faculty of Mathematics, Natural Sciences and Information Technologies, University of

Primorska, SI-6000 Koper, Slovenia ┴

Malaysian Institute of Pharmaceuticals and Nutraceuticals, Ministry of Science, Technology

and Innovation, Halaman Bukit Gambir, 11900 Bayan Lepas, Pulau Pinang, Malaysia

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 55

ABSTRACT

Increased reports of oseltamivir (OTV) resistant strains of the influenza virus, such as H274Y mutation on its neuraminidase (NA), have created some cause for concern. Many studies were conducted in the attempt to uncover the mechanism of OTV resistance in H274Y NA. However, most of the reported studies on H274Y only focused on the drug-bound system, but direct effects of the mutation toward NA itself prior to drug binding still remain unclear. Therefore, molecular dynamics simulations of NA in apo form, followed by principal component analysis and interaction energy calculation, were performed to investigate the structural changes of the NA binding site as a result of H274Y mutation. It was observed that the disruption of the NA binding site due to H274Y was initiated by the repulsive effect of Y274 on the 250-loop, which in turn altered the hydrogen bond network around residue 274. The rotated W295 side chain caused the upward movement of the 340-loop. Consequently, sliding box docking results suggested that the binding pathway of OTV was compromised due to the disruption of this binding site. This study also highlighted the importance of the functional group at position C6 of sialic acid mimicry. It is hoped that these results could improve the understanding of OTV resistances and shed some light on the design of a novel anti-influenza drug.


2

Page 3 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION Influenza, also referred to as the “flu”, is a disease whose transmission could reach pandemic levels. The worst influenza pandemic was estimated to have killed more than 50 million people in 19181. Influenza is caused by the influenza virus, which has evolved and mutated over time, posing serious challenges in developing effective therapies. In the last century, at least three influenza pandemics have been reported: the Spanish Flu (H1N1) in 1918, the Asian Flu (H2N2) in 1957 and 1968, and the Hong Kong Flu (H3N2). In 2009, a pandemic swine flu (H1N1) was reported to have killed approximately 18,500 people. Compounding this serious health issue is the threat from the avian influenza virus, (H5N1), which is originally carried among birds and produces a high fatality rate when transmitted to humans. In 2011, there were 34 deaths out of 62 cases of H5N1 reported from Bangladesh, Cambodia, China, Egypt, and Indonesia2. Most recently in 2013, a novel avian influenza virus, H7N9, which was not previously known to infect humans, resulted in 37 deaths out of 132 patients in China3-5. The glycoproteins on the surface of the influenza virus, i.e. hemagglutinin (HA) and neuraminidase (NA), have become targets for the development of anti-influenza therapies, as reviewed by Itzstein et al.6. Up to 2013, there have been 18 HA and 11 NA subtypes of influenza7. HA is the main target for vaccine development due to its high immunogenicity, although some small molecule HA inhibitors have also been reported. The major challenges in targeting HA are due to the shallow, less energetic, and accumulated mutations at the receptor binding site8, 9, 10. But the structural conservation of the NA binding site, the important role of NA in viral infections, and the availability of NA crystal structures in the Protein Data Bank (PDB) provide some advantages for NA as the most targeted protein for inhibitor development in


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 55

anti-influenza drug design6. Tamiflu (a pro-drug for oseltamivir carboxylate, OTV) and Relenza (zanamivir, ZMR) are two examples of the recently FDA approved NA inhibitors. Although Tamiflu is the world’s largest stockpile of the anti-influenza drug and most commonly used for influenza treatment, the increased reports of OTV resistant strains (e.g. H274Y, R292K, N294S, and the recent H7N9 China strain) have created cause for concern11-16. Thus, there are considerable efforts in the pursuit of antiviral small molecules as OTV alternative via ligand and structure based drug design approaches, natural products and analogue synthesis17-19. Developing a deep understanding of OTV resistance mechanisms, especially at the molecular/atomic level, has become crucial. There have been many studies (mutagenesis20, crystallography21, and molecular dynamics simulations22-27) conducted in an attempt to uncover the mechanism of OTV resistances due to H274Y mutation (see Figure 1). In a mutagenesis study20, it was observed that residues with a bulkier side chain (i.e. tyrosine, phenylalanine) reduced the NA sensitivity toward OTV by posing steric hindrance for reorientation of the nearby E276. The rearrangement of E276 was required to form a proper hydrophobic pocket for OTV6, 28. However, this effect was not observed for smaller residues or ligands with glycerol moieties attached to them, such as ZMR and 4-amino-DANA. Interestingly, OTV affinity in the H274Y mutant of recombinant H3N2 was not affected20, which suggests that the differences in residues surrounding position 274 in the N2 and N1 subtypes play an important role in the resistance to OTV. The results obtained from the mutagenesis study20 were supported by the NA crystal structures of the H274Y mutant from the H5N1 virus in complex with OTV and ZMR21. As observed in the crystal structure of the H274Y-OTV complex (PDB ID: 3CL0), the bulky side chain of Y274


4

Page 5 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


pressed the carboxylate side chain of E276 ~2 Å deeper into the binding site, as compared to that in the wild type (WT). This resulted in the disruption of the hydrophobic pocket, which could not accommodate the pentyloxy group of OTV. However, the disruption did not affect the two ZMR hydroxyl groups’ ability to maintain the interaction with the carboxylate group of E276. This observation agreed with the molecular dynamics (MD) studies of the H274Y-OTV complex system, where the conformational change of E27626, 29 and the loss of the bidentate salt bridge between E276 and R22426 were also observed. The MD results showed that Y274 forced the rotation of the carboxylate group side chain of E276 by 115°, thus altering the hydrophobic subpocket. Due to the altered binding site, the pentyloxy group of OTV was rotated by 125°, reducing the binding energy by ~5 kcal/mol23. In another MD simulation study of the OTVH274Y complex system, it was proposed that the origin of OTV resistance was due to the infiltration of water molecules into the binding site27. MD studies of various mutants (e.g. H274Y, N294S, and R292K) in complex with OTV22,

30

indicated that each mutation has resulted in different conformational changes at the active site. Therefore, the mechanism of drug resistance in various mutants should not be generalized. Furthermore, the differences between the buried surface areas of OTV’s pentyloxy in the WT and in H274Y were not in agreement with those suggested from the crystallographic data21. The fact that the active site of NA is composed of many charged residues—i.e., arginine triads (R118, R292, and R371), E118, E119, E276, E277, D151, R152, and R224—the electrostatic potential has been suggested as the driving force of ligand diffusion and its stability in the active site31-33. In fact, H274 did not directly interact with the ligand, despite being one of the framework residues to support the structure of the catalytic site34, 35. Therefore, studies that only focused on the endpoint of ligand protein interaction to explain the effects of this mutation might


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 55

not be sufficient. Normal mode analysis and MD simulations of NA with and without the presence of OTV to investigate the binding kinetics demonstrated that the binding site in the H274Y mutant was more rigid than that in the WT25. Thus, OTV could exhibit stronger binding to WT than to mutant type (MT). Unfortunately, the reported timescale (1 ns) might not be enough to suggest a strong conclusion. Although the disruption of OTV binding in the H274Y NA active site is extensively investigated, Le et al.24 have highlighted the need to further investigate the effect of H274Y in the process of drug binding based on the observation of the binding funnel disruption. The reported binding kinetics parameters21 also indicated that OTV and ZMR were slowly associated with NA in MT compared to WT. The changes of drug-receptor interactions between OTV and NA binding sites due to H274Y have been studied20-23, 26, 27, 29, 30. However, the direct effect of H274Y mutations toward the dynamics behavior of the NA binding site itself, albeit prior to the drug binding, still remains unclear. To the best of our knowledge to date, most of MD simulation studies on the resistance mechanism of OTV have dealt only with ligand-bound (holo) systems22, 23, 26, 27, 29, 30, 36, 37, with a few exceptions24, 25. Here we report the results from MD simulations of various NA systems in the apo form. Our goal was to investigate the specific role of H274 in the WT. In addition, we were interested in the direct effects of the substitution of histidine with tyrosine at this position on the integrity of the NA binding site and the correlation of the effects toward OTV binding. METHODS System Preparation


6

Page 7 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Six apo-NA MD systems were prepared (see Table 1). NA structures used in the simulations were taken from PDB (www.rscb.org38). All cysteine residues in NA formed disulfide bonds. Hence, the notation of CYS for cysteine was changed into CYX according to the amino acid type in AMBER. Native wild type N1 systems, based on PDB ID 2HTY39, were built for each protonation state of histidine, namely, HID (protonated δN), HIE (protonated εN), and HIP (both of δN and εN protonated). The native H274Y N1 mutant system labeled MT in Table 1 was built by removing the co-crystallized ligand ZMR from the PDB structure 3CKZ21. In order to ascertain that the conformational changes in the MT system were induced by the H274Y mutation and did not occur from the removal of the ligand from the co-crystal structure, we carried out another two simulations using different initial structures. The simulations were: (1) built from the apo form of N1 from the crystal structure of the WT N1 (PDB ID: 2HTY) and then mutated H274 to tyrosine (WT→MT in Table 1), and (2) built from the holo form of H274Y N1 (PDB ID: 3CKZ) and then mutated back to the WT (i.e., Y274 into histidine) (see HY→WT in Table 1). The mutation of amino acid residues in WT→MT and HY→WT was accomplished using the Build and Edit Protein module in Accelrys Discovery Studio40. Conserved calcium ion and crystal water molecules were included in the simulations. A box of TIP3P water was added to solvate each system using the tleap program resulting in a rectangular box with a size of ~80 Å3. The minimum distance between the protein and the edge of the box was set to 10 Å. The system was subsequently neutralized and introduced to a 150 mM NaCl salt bath. Molecular Dynamics Setup


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 55

Minimization and MD simulations were performed using the PMEMD module in AMBER1141 with 32 core CPU and the AMBER FF03.r1 force field42. Each system was first subjected to 1,000 steps of steepest descent minimization. Then 2,000 steps of conjugate gradient minimization with 5 kcal/molÅ2 of harmonic restraints were applied. The restraints were gradually released from all heavy atoms, including on the side chains, to only backbone atoms with the same minimization protocol. A final 5,000 steps of unrestrained conjugate gradient minimization were performed in order to remove any steric clashes. Following that, the system was gradually heated to 310 K (0-100 K; 100-200 K; 200-310 K for 50 ps each) in the NVT ensemble using harmonic restraints on the backbone atoms. Next, 1 ns of NPT equilibration was performed, where the harmonic restraints on the backbone were gradually decreased by 1 kcal/molÅ2 per 250 ps. Then a 20 ns production run in the NPT ensemble was performed with all hydrogen atoms constrained using the SHAKE algorithm43. The temperature was controlled with a Langevin thermostat44 with a collision frequency of 1 ps-1, and the pressure was controlled using a Berendsen barostat45 with a coupling constant of 1 ps and a target pressure of 1 bar. The time step used during the production phase was set to 2 fs. A non-bonded cut-off value of 10 Å was used, and the long range electrostatics were treated using Particle Mesh Ewald46. Postsimulation analysis was conducted using Ptraj module in AmberTools 1.541. Free Energy Calculation Pairwise interaction energy was calculated using the MM/GBSA scheme in MMPBSA.py version 1347. In the MM/GBSA calculation, snapshots were taken per 1 ns from the 20 ns of production trajectory. This calculation compared the free energies of multiple conformations, from continuous snapshots, to determine their relative stability. The free energy associated with conformational change, from state A to state B, was calculated using (1):

∆GA→Bsolvated = ∆GB,solvated - ∆GA,solvated

(1)


8

Page 9 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The change of free energy of each term on the right-hand side of Equation (1) was estimated using the following equation: (2)

∆ = + ∆ − !"

(3)

∆ = ∆#$ + ∆%&

The gas-phase energies (Egas) were calculated from the molecular mechanics (MM) energies, i.e. van der Waals and electrostatic interactions. ∆Gsolvation composed of polar and non-polar solvation which were calculated using implicit solvent model (Generalized Born, GB) and solvent-accessible surface area (SASA). S is the entropic contribution. The energies was estimated according to the averages of multiple conformations taken from MD trajectory: ∆ ≅ 〈 〉 + 〈∆ 〉 − 〈!" 〉 *

*

(4)

-

+ + = + ∑+ .* , + + ∑.* ∆, − + ∑.* !,"

where i is the particular snapshot and N is the total number of snapshot included in the calculation. Pairwise decomposition calculates the interaction energy between pairs of residues in the system.

Principal Component Analysis


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 55

The 20 ns MD trajectory of each WT, MT, WT→MT, and HY→WT were superimposed to the first frame using ptraj module in AmberTools 1.5. Principal component analysis was performed on each CA atom of the protein backbone from these MD snapshots using ProDy48 in VMD software49. Sliding Box Docking Snapshots of the WT and MT systems were extracted for every 2 ns of the production run, giving a total of 11 snapshots. Each of these snapshots was reoriented so that its active site was directed toward the x-axis. An initial searching space (grid box) was created to cover the top surface of NA with the size of 10 x 48 x 48 Å in x, y, and z, respectively (adapted from Tran et al.50). The sliding grid box from the top surface to the final active site was accomplished by modifying the x-coordinate of the grid box for every -2 Å, and it yielded six total boxes, denoted Box-1 through Box-6 (see Figure S1). OTV and ZMR were docked to boxes 1 through 6 using AutoDock Vina51 with the number of docking outputs set to 20. Docking poses from Box-1 and Box-2 were merged to represent the ligand’s conformations at the entrance of the binding pathway. Ligand’s conformations at the intermediate and the final active site regions were represented, respectively, by merging Box-3 with Box-4 and Box-5 with Box-6. The results were visualized using Accelrys Discovery Studio Visualizer52. The plots of the docking results of each snapshot are shown in the Supporting Information section (Figures S4 and S5). The Connolly molecular surfaces were applied to each receptor snapshot using probes with a radius of 1.4 Å.

RESULTS AND DISCUSSION


10

Page 11 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The role of native-histidine at position 274 in WT In a mutagenesis study, Wang et al.20 demonstrated that any substitutions at position 274 might change the affinity of OTV, ZMR, and 4-amino-DANA toward NA. This implied that H274 itself might have a distinct role in the NA binding site. It is noted that under physiological pH, histidine’s side chain has two nitrogen atoms that can be protonated, δN and εN. From the crystal structure of N1 WT in complex with OTV (PDB ID 2HU439), both δN and εN of H274’s side chain have possible interactions at a distance of ~3 Å, i.e. with P245, G248, N294, and W295 (See Figure 2). Although the protonation state of histidine can be predicted based on pKa calculation and the optimum hydrogen bond pattern in crystal structures, MD simulations of a solvated system could offer a more conclusive result53. In order to investigate the role of native histidine in all possible states on the structure and dynamics of WT NA in apo form, we simulated these three protonation states in WT systems: 1) all δN histidines were protonated (HID), 2) εN were protonated (HIE), and 3) both of δN and εN were protonated (HIP). The RMS deviations (RMSD) and RMS fluctuations (RMSF) calculated for the protein backbone atoms over the 20 ns of production trajectory are presented in Figures 3 and 4, respectively. The RMSD plot (Figure 3) indicates that the optimum overall stability of the NA structure was achieved in the HIE system, the lowest RMSD among the WT systems. Figure 4 shows that the loop regions, namely 150, 250, 270, 300, 340, and 380, demonstrated high RMSF values. Higher fluctuations for the N and C termini were due to the fact that these regions were not restrained. The high fluctuation of the 150-loop was expected since it was recognized previously as a flexible loop in Group-1 NA54, 55. Among the three WT systems, this loop showed the lowest fluctuations in HIE. In addition, the HIE system also generally exhibited the lowest fluctuations for loops around H274, specifically the 250, 270, 300 and 340 loops, suggesting that the HIE system has the most


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 55

optimum interactions among these loops in the WT. Nevertheless, the participation of the amino acid side chains on the dynamics of these loops was not represented by the RMSF values above. Therefore, the 20 ns average structure of each system was generated and compared. Figure 5 shows that the protonated εN of H274 in the HIE and HIP systems acted as a hydrogen bond (Hbond) donor to the residues of the 250-loop, i.e., with P245 and G248 (Figures 5b and 5c). These H-bonds did not occur in the HID system since εN of H274 was not protonated. Instead, its protonated δN acted as an H-bond donor to the backbone of R292 on the 300-loop (Figure 5a). An additional H-bond in the HIE system was formed between δN of H274 as the acceptor and the backbone of Y275. Although an average structure may reflect the dynamics of the MD system, it does not show the time-dependent structural changes throughout the simulation. For this reason, time evolution snapshots over 20 ns of each production trajectory were captured every 1 ns. The snapshots focused on the loop regions around residue 274 (Figure 5e-g). In the HID system, the high motion of the 250-loop, as shown in Figure 5e, corresponded well with the RSMF values. As discussed before, the 250-loop in the HID system did not form an H-bond with the side chain of H274, thus it became destabilized and deviated from its initial position (Figures 5a and 5e). In the HIE system, there was no significant divergence among the snapshots throughout the 20 ns (Figure 5f). This agreed with the low profiles of the RMSD and RMSF values. However, in the HIP system, although its 250-loop showed better stability compared to the HID, noticeable movements of the side chains of H274 and R292 were detected over the 20 ns simulation (Figure 5g). Therefore, the RMSD and RMSF plots (Figures 3 and 4, respectively) and the structural dynamics shown in Figure 5 have suggested that the most favorable protonation state of H274


12

Page 13 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


was HIE, when the εN of the imidazole nitrogen was protonated. Moreover, it was noted that there were seven histidines in the structure of NA: two of which were located around the H274Y mutation point (H274 and H296), the other two around the 150-loop (H144 and H155), and the rest located farther from the binding site (i.e. H126, H184, and H412A). The pKa of H126, H144, H155, H184, H274, H296, and H412A calculated using PDB2PQR56 server, were 5.97, 6.32, 6.20, 4.37, 5.19, 6.16, and 4.81, respectively57, 58. Thus at pH 7, these seven histidines were predicted to be in neutral state, either HID or HIE. However, the final protonation state of histidine was determined based on the lowest RMSD achieved in the system53. Based on the RMSF plot, the stable profile of the 150, 250, 270, 300, 340, and 380 loops in the HIE system might be due to the occupations of H144, H155, H274, and H296 in the protonated εN state. In addition, the residual fluctuations of H126, H184, and H412A in the HIE system were comparatively lower than in the other systems (Figure 4). For this reason, it is suggested that the most preferable state for all histidines in the NA structure is when only the εN is protonated. Subsequently, we simulated H274Y mutant type NA (denoted as MT) for comparison with all the wild type systems above. The MT system was built based on the crystal structure 3CKZ, which has the best resolution (1.9 Å) among all the reported crystal structures of N1-H274Y in the PDB. Since 3CKZ is a holo form of NA in complex with ZMR, ZMR was removed from the NA structure prior to simulation. All histidines were assigned as protonated εN, because it was demonstrated earlier that this state was the most preferable in the WT simulated systems. In general, the RMSD values of the protein backbone in the MT system over 20 ns were higher than those in the HIE system (Figure 3). The residual RMSF plot of the MT system indicated that the loops around the H274Y mutation point, i.e., the 250, 300, and 340 loops, fluctuated more than those in the HIE system (Figure 4). The 20 ns average structure of the MT system showed


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 55

that no H-bond formed between the side chain of Y274 and its surrounding loops (Figure 5d). The motions of the 250-loop (Figure 5h) corresponded well with the RMSF values. However, it is predicted that the roots of the 250-loop motions in the HID and MT systems are different. In the MT system, tyrosine, being a larger moiety compared to histidine, appeared to compel the 250-loop away from the 270-loop (Figures 5d and 5h). The steric repulsion between the phenolic oxygen of Y274 and the carbonyl oxygens of both P245 and S246 could be why the 250-loop moved apart from the 270-loop. The results also indicate that even a small change in the atomic composition in one amino acid residue, e.g., one hydrogen atom difference in the various protonation states of H274 or the mutation of H274 into Y274, has the ability to disrupt the stability of its surrounding residues. It is also interesting to note that from the last snapshot of the MT system, the side chain of R292 was leaning toward E276 (Figures 5d and 5h). The time-dependent mechanism of the E276-R292 formation will be further discussed below. Effect of H274Y Mutation Conformational changes in the vicinity surrounding position 274 In order to ascertain that the conformational changes in the MT system were induced by H274Y mutation and not from the ligand’s removal from the co-crystal structure, MD trajectories of the WT and MT systems were compared to those of the manually-mutated systems, WT→MT and HY→WT (see Methods and Table 1). Figures 6 and 7 show that, in general, systems built from the holo form crystal structure (i.e., with the PDB structure 3CKZ; MT and HY→WT) had higher conformational changes from their


14

Page 15 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


backbone motions than those based on PDB ID 2HTY (WT and WT→MT). This result was expected, since 3CKZ is the ligand-bound NA structure. Removing the ligand caused these systems to find a new equilibrium state. Higher RMSF values were noted in regions representing the 150, 250, 270, 300, and 340 loops for the H274Y systems compared to the WT systems (Figure 7). As for HY→WT system, despite being a wild type, the fluctuation in the 150-loop region might be due to the removal of ZMR from the holo form as previously discussed. In the holo form 3CKZ crystal structure, the 150-loop was in closed conformation to stabilize its ligand, ZMR59. This loop might have transformed from closed to open form due to the absence of its ligand, which resulted in higher RMSF in this region for MT and HY→WT. It is worth noting that the higher RMSF of the 340-loop region was only observed in the H274Y systems (MT and WT→MT) and not in the wild types (WT and HY→WT). To the best of our knowledge, a high fluctuation of the 340-loop due to H274Y has not been reported so far. Superimposition of the time evolution snapshots at the interval of 4 ns throughout the 20 ns of MD simulations (Figure 8) also showed that significant conformational changes occurred at the 250, 300, and 340 loops (Figures 8b and 8c), which corresponded well with the RMSD and RMSF values. The amino acid compositions for these loops are as follows: 250-loop (residues 244-250), 270-loop (residues 268-275), 300-loop (residues 293-300), 340-loop (residues 328348), and 150-loop (residues 147-152)60. In all the H274Y systems, the 340-loop was disturbed, especially in the MT system; i.e., this loop was found to shift upward (Figures 8b and 8c). Similarly, the side chain tilting of R292 toward E276, of N294 toward the 250-loop, and the rotation of the W295 side chain were observed in all H274Y systems. Specifically in the MT system, Y274 was found to flip back toward N272 at 20 ns (Figure 8). In contrast, there were no


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 55

major conformational changes observed in the WT systems except at the C- terminal as evidenced from the superimposition of the snapshots (Figure 8a and 8d). Hydrogen bond network in the vicinity surrounding position 274 H-bond formation in the vicinity surrounding residue 274 (250, 270, and 300 loops (Figure S2) was monitored throughout the 20 ns of simulation using the distance and angle cutoff of 3.00 Å and 120.00°, respectively. All detected H-bonds (>1% of occupancy) are listed in Table S1. Hbonds that showed significant divergence due to H274Y mutation are listed in Table 2. It is interesting to note that the 250, 270, and 300 loops in the WT possessed an extensive hydrogen bonding network. H-bonds that were involved in the interactions between the 250 and 270 loops (T242-Y275, S246-E276, and Y252-Y273), between the 250 and 300 loops (N247-N294), and between the 270 and 300 loops (H272-H296) were found to be less persistent in all H274Y systems than in the wild types. High occupancy of the H-bond between T242 (located at the βstrand immediately before the 250-loop) and Y275 (member of 270-loop) in the wild types (WT and HY→WT systems) resulted in the close proximity between the 250 and 270 loops. This Hbond has 64.8% occupancy in the WT, but decreased to 51.4% when the WT was mutated into H274Y, WT→MT. On the other hand, its occupancy in the MT system was 18.6%, but it increased to 65.0% when the MT was mutated back into the wild type HY→WT. The decrease in the interaction between T242 and Y275 in the H274Y systems resulted in the 250-loop to begin losing the interaction with the 270-loop, thus moving away from the 270 and 300 loops. A reduction in H-bond occupancy due to H274Y mutation was also observed for the H-bonds between R292 and N294 in all H274Y systems. On the contrary, high occupancy of the H-bond between E276 and N294 was found in the MT system. The H-bond of R292-D293 was formed in


16

Page 17 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


all H274Y systems (at 2 ns in MT and at 15 ns WT→MT systems) prior to the formation of the H-bond of E276-R292 (at 9 ns and 18 ns in MT and WT→MT systems, respectively). The time evolution of the H-bond formation between R292-D293-N294-E276 pointed to the fact that the diminishing ability of N294 to form a H-bond with R292 in H274Y systems has instead led to the formation of a H-bond between E276 and R292, which in turn, was assisted by the hydrogen bonding between R292 and D293 H-bonds (Table 2). Therefore, from the present results, it is clear that N294 also plays an important role in stabilizing R292 to maintain the integrity of the NA binding site. Table 2 shows that the H-bond formation between the side chains of E276-R292 distinctively occurred in all H274Y systems, including the MT and WT→MT systems, but not in the wild type systems. As for the HY→WT system, E276-R292’s H-bond was observed at a very minimum occupancy. In the MT system, this H-bond was formed since the first nanosecond despite its low occupancy. It then reached an occupancy of greater than 50% after 9 ns. As for the WT→MT system, it was only formed after 18 ns. This result reflected that H274Y mutation caused the H-bond of E276-R292 to form faster in its native crystal structure (3CKZ) than in the WT→MT system. It is noted that the MD simulation of the WT→MT system was started from WT conformation (apo form 2HTY). The newly introduced Y274 in the WT→MT system caused its WT conformation to adapt slowly until it finally resulted in the H-bond formation between the side chains of E276-R292. The positively charged arginine has the ability to interact with the negatively charged glutamic acid through strong electrostatic interactions, e.g., H-bond and salt bridge. With the exception of the work by Rippol et al.30, where a salt bridge formation between R292 and E276 was found in the H274Y-OTV complex system, no other publications


17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

on MD studies of H274Y of N1 subtype22-27,

61

Page 18 of 55

have reported or highlighted the formation of

E276-R292 salt bridge. The movement of the side chain of R292 toward E276 (Figure 9) in H274Y systems could contribute to the OTV resistance in NA, especially in the N1 subtype. R292 is located adjacent to the hydrophobic pocket of OTV formed by E276 and R224 (Figure 9a). OTV possesses a slow binding characteristic because of the hydrophobic subpocket6,28 formation requirement. Therefore, the formation of an H-bond between E276 and R292 created a polar environment, which might be unfavorable to the binding of the OTV pentyloxy side chain into the subpocket (Figure 9b). Inspection of the complex ZMR-H274Y and OTV-H274Y crystal structures (3CKZ and 3CL0, respectively21) revealed that there was no H-bond formed between the side chains of E276-R292. However, the kinetic studies of OTV have indicated that it has a lower association constant and a higher dissociation constant in the H274Y mutant than in the wild type21. Therefore, it can be predicted that the formation of the H-bond of E276-R292 prior to ligand binding could result in a longer time needed for the formation of the N1 hydrophobic subpocket in accommodating OTV. Pairwise decomposition of interaction energies between E276-R292, R224-P245, T242-Y275, W295-H296, N294-R292, N294-N247, R292-D293, and E276-R224 (the important H-bonds formed in the simulations), were calculated using the MM/GBSA method (Figure 10 and Table S2). Weaker energies of R224-P245, T242-Y275, and N294-N247 occurred in all H274Y systems than in the wild types corresponded well with the increased distance between P245 (part of the 250-loop) and R224; between T242 and Y275 (of the 250 and 270 loops, respectively); and between N247 and N294 (of the 250 and 300 loops, respectively). This could possibly be


18

Page 19 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


due to the flexible motion of the 250-loop as a result of Y274 steric repulsion observed in the H274Y systems. Based on the timeline of H-bond occurrences in the H274Y systems (Table 2), the interaction between the side chains of E276-R292 was initiated by decreasing the H-bond occupancy between N247-N294. As a result, N294 lost its interaction with R292. Thus, the destabilizedguanidine side chain of R292 finally approached E276 with the support of the H-bond of R292D293. These observations were in agreement with their energy profiles. The weakened energy of N294-R292 was observed in the MT and WT→MT systems. In contrast, stronger interactions of R292-D293 and E276-R292 were observed in all H274Y systems compared to the interaction in the wild types. However, compared to the MT, the weaker energy of E276-R292 and R292D293 in the WT→MT were expected since these H-bonds formed faster in the MT system. In addition, relatively weaker interaction energy between R224-E276 was observed in the H274Y systems, i.e., MT and WT→MT systems, than in the wild types (Table S2). This result agreed with the previous observation regarding the disruption of the salt bridge between R224E27612, 26, which was attributed to H274Y mutation. Although the interaction was quite weak, its relatively high H-bond occupancy in both WT and mutants showed that it did not completely dissipate. This also corresponded with the consistent H-bond of R224-E276 in the H274Y system23 that was previously observed. From the H-bond analysis and pairwise energy decomposition, it is clear that the repulsive steric effect of Y274 toward the 250-loop led to the disruption of many H-bond networks surrounding residue 274. Such disruption could be observed between the 250 and 270 loops (T242-Y275,


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 55

S246-E276, Y252-Y273); between the 250 and 300 loops (N274-N294); and between the 270 and 300 loops (N272-H296). Dihedral analysis In the crystal structure of WT in apo form (PDB ID 2HTY39), W295 appeared to be stabilized by H296 through the π-π stacking interaction61. However, the present MD simulation showed that at 20 ns the side chain of W295 in H274Y systems was rotated, and so it lost its interaction with H296 (Figures 8b and 8c). Further dihedral analysis showed that the side chains of W295 and H296 generally rotated more frequently in the H274Y systems than in the wild types (Figures 11a, 11d and 11e), leading to the disruption of the H-bond between Y272-H296 (Table 2). This in turn destabilized the side chain of W295, as indicated by the gradual changes of its dihedral angle after 7 ns (Figures 11a and 11d). In addition, the interaction energies between W295-H296 also decreased in the H274Y systems (Figure 10 and Table S2). This concurred with the rotation of the W295 side chain in all the H274Y systems, which appeared to disrupt the π-π stacking between W295-H296. It is also noted that dihedral angle analysis on the side chain of H/Y274 showed that the side chain of Y274 in the MT system appeared to flip back toward N272 at 20 ns (Figures 8c, 11c). The rotation of Y274 side chain reached a new stable state at 18 ns, which corresponded to when the H-bond formation between N272-Y274 occurred (Table 2). Interestingly, this incident also coincided with the flipping of H296 in the MT system at 18 ns (Figure 11a). It is worth noting that the inclination of Y274 side chain to reorient toward N272 was also observed in the crystal structure of H274Y of the N3 subtype (PDB ID: 4HZY)61.


20

Page 21 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Essential dynamics of N1 binding site From this study, H274Y has been observed to affect the structural dynamics around the mutation point. High fluctuations of loops surrounding the H274Y-NA binding site, i.e., the 150, 250, 270, 300, and 340 loops, were observed throughout the 20 ns simulation. Interestingly, when the H274Y crystal structure was mutated back into the wild type (HY→WT), the RMSF of these binding site loops showed similar profiles as those in the WT system—except for the 150-loop. The essential movement of the NA backbone was then calculated using principal component analysis (PCA) method48. The motions of the protein backbone and the magnitude of the atomic fluctuations are represented by the directions and length of the arrows (Figure 12). In the WT system, most of the loops at the binding site were stable (Figure 12a). On the other hand, in the H274Y systems, the essential dynamics analysis indicated that the movement of one loop could affect the motions of its neighboring loops. The relatively long arrows on Cα of the 250, 270, 300 and 340 loops in all the H274Y systems indicated high mobility of these regions. Figures 12b and 12c show that the whole 250-loop in the H274Y systems moved away from the 270-loop. This movement was probably due to the steric effect of the bulky side chain of Y274, as suggested earlier. Unlike the 250-loop, not all residues in the 300-loop moved in the same direction. N294 moved toward the 250-loop, while W295 and H296 shifted toward the 340-loop. The arrows indicated upward movement of the 340-loop in all H274Y systems, which corresponded with the final conformations of the 340-loop in the MT system (Figure 8c). Superimposition of the 340-loop between the MT and WT systems at 20 ns also indicated a significant deviation in the motion of the MT system (Figure S3). In the MT system, W295 appeared to be stabilized in a new conformation where its side chain pointed upward and interacted with N344 and G345 of the 340-loop. The attractive force between the pyrrole ring of


21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 55

W295 and G345 was contributed by the π-donor H-bond62. Also, the interaction between the sixmembered indole ring of W295 and N344 (amide-π stacked interaction63) resulted in the 340loop flipping up from its initial position. The consequence of the H274Y mutation toward the 150-loop could be explained in the WT→MT system. In WT system, the 150-loop remained stable in open-conformation, a native state of NA in apo form, with a low magnitude of motions (Figure 12a). An exception of 150loop’s motion in HY→WT was due to the removal of ligand from its closed-conformation (holo form; Figure 12d). However, when H274Y was being introduced to the WT structure i.e. WT→MT, high motion of the 150-loop was observed (Figure 12c), that could be due to the clockwise movement of the 250-loop, as suggested from the magnitude and direction of the arrows in Figure 12c (on the left). This observation indicated that even one mutation such as H274Y has a cascading effect toward the larger region of loops at the active site (starting from the 250-loop to 150-loop). The important role of the loops in NA binding sites to ligand binding/unbinding events has been discussed at length in previous studies24, 36. Therefore, the results of the present study could provide additional insight into slower OTV binding and faster OTV unbinding in the H274Y mutant compared to the wild type21. Recently, long-scale MD simulations of H274Y in complex with OTV36 demonstrated faster opening of the 150-loop in the mutant, which in turn caused OTV to dissociate from the binding site of NA. The essential dynamics of MD systems based on NA in holo form (3CKZ), (i.e., MT and HY→WT), showed high motions of the 150-loop despite HY→WT being a wild type. As discussed before, the 150-loop in the holo form of the N1 crystal structure (e.g., PDB ID 2HU4, 3CKZ, 3CL2, 3CL0) were in a closed conformation to stabilize the bound ligand21,

39, 54

. A


22

Page 23 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


closed form of 150-loop in complex with inhibitors also appeared to be energetically favorable64. Since the ligand was removed in our simulations, the closed conformation of the 150-loop in the MT and HY→WT systems was destabilized and transformed to an open conformation54. Sliding box docking of OTV and ZMR The region around residue 274 has been proposed as the binding funnel for OTV24, 50 and the location where the pentyloxy group of OTV bound to the hydrophobic pocket65. Thus, any changes that occurred at this region would affect the binding of OTV24. The present MD study showed that the presence of Y274 in all mutant systems resulted in the disruption of its surrounding loops, starting from the repulsion of the 250-loop to the upward movement of the 340-loop. Nevertheless, the correlation between these behaviors and the ligand (drug) binding is still unclear. In fact, the mutation of H274Y has been reported to reduce the NA sensitivity toward OTV, but not ZMR16, 20, 21. We suspected that the change of NA active site’s behavior in mutant, would also affect the drug binding process, besides the endpoint interaction of ligand-protein. For this reason, sliding box docking of OTV to the NA binding site was performed and compared to that of ZMR. This docking concept was recently introduced by Andrelly Martins-José66 and has been used to investigate the binding pathway of OTV on the holo form-NA of H5N1 virus50. This method could be useful to predict the possible binding pathway(s) of ligand, from solvent (surface) to the final active site. Despite the rough approximation of ligand binding process, this approach is avoiding initial biases of specifying the starting point and the direction which the ligand follows50. Sliding box dockings of OTV and ZMR were performed on NA snapshots from 20 ns of MD production trajectory of the WT and MT systems (Figures S4 and S5). Docking results of the


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 55

initial snapshot of the WT system (at 0 ns), suggested that both OTV and ZMR were able to enter the active site from two directions: the 340-loop and 250-loop regions (Figures 13a and 13b). The pathways from which the ligand entered (the 340-loop and 250-loop) are referred to as Path-1 and Path-2, respectively. Comparison between the two Connolly molecular surfaces of the WT system, at the initial and final snapshots, showed no notable changes on Path-1 and Path-2 (see inset of WT in Figure 13). These observations reflected the stability of the NA binding site of the WT throughout the 20 ns simulation. A relatively slight change was only observed in the 150-loop region: there was a formation of an extra cavity that provided an additional pathway to the binding site of the WT system as demonstrated by the snapshot at 20 ns (Figures 13c and 13d). The energetic profiles of OTV binding through Path-1 and Path-2 in the WT system were evaluated by plotting the docking scores at each point of the pathways as represented by OTV centroids (Figure S6, left). The results indicated that the binding of OTV through these two pathways was energetically favorable. The docked conformations of these pathways were clustered based on their positions: entrance (the 340 and 250 loop regions in Path-1 and Path-2, respectively; ligand colored in red); intermediate (the 430 and 250 loop regions in Path-1 and Path-2, respectively; ligand colored in orange); and final active site (ligand colored in green). The distributions and the average of docking scores in Path-1 are as follows: the final active site (-4.9 kcal/mol) < intermediate (-4.8 kcal/mol) < entrance (-3.9 kcal/mol) (Figure S6, right). Moreover, OTV binding through Path-2 was also energetically favorable with a similar profile to Path-1 (see Figure S6). This result agreed with that of the sliding box docking of OTV on the holo form of the N1 crystal structure50 where Path-1 and Path-2 were referred to as the tunneling and climbing pathways, respectively50.


24

Page 25 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


It is worth noting that OTV binding on the final snapshot in the H274Y system appeared to be compromised. None of the OTV docked conformations were found at the entrance of Path-1 or Path-2 (Figure 13e). Although none of the ZMR docked conformations were found at the entrance of Path-2, ZMR was docked favorably at the entrance of Path-1 (the 340-loop region) (Figure 13f). The failure for both OTV and ZMR to enter the binding site through Path-2 was probably due to the high motions of the 250-loop in the H274Y systems. However, the ability of ZMR to enter the active site through Path-1 might be correlated with the unaffected affinity of ZMR in the H274Y mutant20, 21. Further investigation into Path-1 in the MT system at 20 ns revealed that the upward movement of the 340-loop created a narrow cavity between the 340-loop and its adjacent 320-loop (see inset of MT surface at 20 ns in Figure 13). The rotation of the Y347’s side chain toward the 250loop resulted in the positively charged R371 being exposed to the surface (Figure 13e). Figure 14 shows docked conformations of ZMR at the entrance of Path-1 in the final MT snapshot. It was observed that the glycerol moiety of ZMR formed at least one H-bond with the exposed R371. It was also noted that the structure of OTV was similar to ZMR, except for different functional groups at the positions of C4, and especially, C6. OTV did not have a glycerol group, which was present in ZMR, located at the C6 atom that formed the H-bond with R371. At the same position, OTV had a hydrophobic pentyloxy group instead of glycerol. Thus, OTV could not form an Hbond with the exposed R371. This was further tested by performing an additional sliding box docking of peramivir (PMV), another NA inhibitor, to the MT snapshot at 20 ns. It was noted that PMV had a similar hydrophobic group with OTV at the C6 position. Interestingly, the results showed that PMV also failed to bind at the entrance of Path-1 (Figure S7). This could strengthen the argument for the importance of the functional group at the C6 position in OTV


25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 55

resistance due to H274Y mutation. Despite the inability to bind at the entrance of Path-1 and Path-2, OTV was still able to reach the active site regions (see the orange and green colored OTV in Figure 13e). This result implied that OTV would eventually reach the active site of the H274Y mutant, as inferred from the lower association rate of OTV binding to the H274Y mutant when compared to its binding to the wild type21. The crystal structure of H274Y in complex with OTV21 revealed the absence of an H-bond between the side chains of E276 and R292. However in our MD simulation, the H-bond formation of E276-R292 in all mutant systems (MT and WT→MT) occurred after 8 ns and 16 ns, respectively (Table 2). Unlike the disruption of the 250 and 340 loops, which were located at the entrance of Path-1 and Path-2, the H-bond of E276-R292 was located at the final active site region. Figure 15 shows that the movement of the R292 side chain toward E276 resulted in an Hbond network between R224-E276-R292. As a consequence, the hydrophobic subpocket intended for OTV became a charged polar subpocket in the H274Y mutant. Furthermore, the best poses of OTV overlaid with the co-crystal structure of OTV in complex with NA (Figure 15a) showed that the hydrophobic pentyloxy group of OTV had an unfavorable interaction with the polar environment at this subpocket. However, ZMR was still able to form a H-bond with E276 and R292 (Figure 15b). Strong electrostatic interactions among R224-E276-R292 in the H274Y mutant would have resulted in a slower formation of a hydrophobic subpocket for OTV binding. This finding might be correlated with the low association rate of OTV in the H274Y mutant21. FURTHER DISCUSSION Most MD studies on the H274Y mutant mainly focused on the changes toward OTV binding at “end-point interaction” with the active site of NA22, 23, 25, 27, 29, 30, 36. However the present study


26

Page 27 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with an MD simulation of NA in apo form aimed at investigating the dynamics of NA prior to ligand binding. The manually mutated structures HY→WT and WT→MT clearly showed the expected overall behavior of the wild and mutant types, thus demonstrating that our simulations indicated that these behavioral changes were because of the H274Y mutation and not a result of removing the ligand from the crystal structures. In the case of OTV resistances, we speculated that H274Y might have an effect on the structural dynamics of NA itself apart from the changes in drug-receptor interactions at the active site. As a result of the H274Y mutation, several distinct motions of loops surrounding the NA binding site were observed: the repulsion of the 250-loop; the upward movement of the 340-loop; and the high motions of the 150-loop contributing to the overall fluctuations of the surrounding loops. In 2012, Woods et al. reported the PCA of 100 ns simulations of NA mutants (the dual mutations of I222R and H274Y) in complex with OTV37. High 150-loop motion was suggested as the major factor to OTV resistances. Although not discussed by the authors, we noted that the PCA results, which were highlighted in their supplementary section, also revealed relatively high motions of the 250, 300, and 340 loops37, which agreed with our finding that the H274Y mutation affected the dynamics of the 250, 340, and 150 loops. Moreover, our proposed binding pathway of OTV and ZMR onto the NA active site in the WT was similar to that which was reported by Tran et al.50, but with the exception of the observed addition of the 150-loop entrance at the final snapshot of the WT systems. Subsequently in 2013, Woods et al.36 extensively studied the OTV unbinding and rebinding events on NA mutants. The role of the 150-loop in releasing OTV from the active site was highlighted. It was noted that the opening of the 150-loop in the H274Y mutant occurred at the same time that OTV moved out of the binding site. In addition, the simulation was able to


27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 55

reproduce previous findings related to the effect of H274Y, specifically relating to the hydrophobic subpocket for OTV decreasing in size and the infiltration of water molecules to force OTV out of the binding pocket27. However, one pertinent question on how the loop motions could help or hinder the ligand binding and its release remains to be answered. Thus, our simulations were purposely carried out in the absence of OTV to understand further what happened before the binding of OTV to NA and how the H274Y mutation affected the dynamics of NA prior to the binding. As a result, the roles of the 250 and 340 loops were highlighted. These two loops were proposed as the entrance of the OTV binding pathways. From the sliding box docking results, the disrupted conformations of the 250 and 340 loops might hinder OTV binding at the entrance of the binding site. This observation might provide an alternative explanation for the lower association rate of OTV into H274Y as compared to the WT21. It is worth noting here that a recent apo form crystal structure of the H274Y mutant, an N3 subtype, also revealed that the 250-loop has a high temperature factor profile61. The importance of 250loop in ligand binding might also be related with an OTV-resistance phenotype that was resulted from the deletions of residues 245-248 (here denoted as 250-loop) in NA recombinant protein67. E276 and R224 are known to form the subpocket for hydrophobic side chains of OTV and PMV68. We observed that the formation of the H-bond between E276 and R292 was not instantaneous, but occurred only after 9 ns in the MT systems and after 18 ns in the WT→MT systems. The formation of these electrostatic interactions (hydrogen bond and salt bridge) might restrain the formation of a proper hydrophobic subpocket, which is required for OTV binding12. For this reason, the low association rate of OTV in the H274Y mutant might also be related to the higher energy barrier in the hydrophobic subpocket formation because of the interaction between E276 and R292. The H-bond of E276-R292 was previously observed in a calcium-


28

Page 29 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


deficient crystal structure of a N9 subtype69. From their MD simulations, Lawrenz et al.70 also observed the strong electrostatic interaction between R292 and E276 side chains in a calciumdeficient NA system. It is noted that the 340-loop is one of the loops present at the calciumbinding site of NA. The absence of a calcium ion could destabilize the motions of the surrounding loops. Hence there might be a correlation between the high motions of the 340-loop and the tendency of the R292 side chain to form an H-bond with E276. Currently, the crystal structure of the H274Y N1 subtype is only available in holo form (PDB ID 3CKZ and 3CL0). Careful observation of the atomic coordinates around the 340-loop in PDB 3CKZ structure revealed an atypically long peptide bond (2.07 Å) between S342 and S344 (S342 and S343 in PDB 2HTY numbering) at the 340-loop region. This peptide bond is longer than the normal partial-double bond of 1.32 Å71. Interestingly, another crystal structure of H274Y in complex with OTV (PDB ID 3CL0, with a lower resolution than 3CKZ) also revealed an unusual length for this peptide bond: 2.24 Å. In addition, the temperature factor of S34421 is comparatively higher than its surrounding residues (see Figure S8). The residue 343 is not evident in PDB 3CKZ and 3CL0. However, this might not be due to the missing residue in the sequence of 3CKZ. Sequence alignment between PDB 3CKZ21 and 2HTY39 confirmed that their sequences are identical, except at position 274. Moreover, a positive contour from the electron density maps of the crystal structures was between a Cα of S342 and S344. This implied the flexibility of the 340-loop in the H274Y crystal structure, which could not be fully resolved during the refinement process. This observation could explain the odd peptide bond of S342S344 in 3CKZ and 3CL0. Nevertheless, in our simulation, high fluctuations of the 340-loop consistently occurred in all H274Y systems; MT and the manually mutated WT→MT. Therefore, we believe that the high motions of the 340-loop was due to H274Y, not the missing


29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 55

peptide bond between S342 and S344 in PDB 3CKZ. This could strengthen our conclusion regarding the mutation effect of H274Y toward the structural dynamics of NA, which could not be explicitly explained in the crystal structure. The timeline mechanism of this observation is shown in Figure S9. During the 20 ns of simulation, it was noted that the magnitude of the upward movement of the 340-loop in WT→HY was not as high as in MT. In the MT snapshot at 20 ns, W295 formed interactions with members of the 340-loop, i.e. N344 and G345 (using PDB 2HTY numbering). We have extended the simulations up to 40 ns and found that the distances between CA atoms of residues N344 and G345 in MT and WT→MT shortened after 40 ns of simulation. The distances for N344 and G345 were of 5.1 Å and 4.5 Å in MT and WT→MT, respectively at 20 ns, then decreased to 4.7 Å and 2.8 Å at 40 ns (Figure S10a). This reflected the closer conformation of the 340-loop in WT→MT than in MT after 40 ns of simulation. In addition, Figure S10e shows that the side chain of W295 in WT→MT formed a H-bond with G345 at 40 ns similar to that in the MT system at 20 ns. This H-bond was not formed in WT→MT at 20 ns (Figure S10d). The occurrence of flip-up conformation of 340-loop in MT system was further examined by extending the simulation up to 30 ns. The time evolution snapshots (Figure S11a) and the RMSD profile (Figure S11b) of 340-loop’s conformation show that the flip-up conformation was stable until 30 ns. This might strengthen the indication of 340-loop’s high motion in the crystal structure of H274Y mutant of N121. In addition to that, a separate simulation of MT system using different random seeds was also performed up to 63 ns (MT-II system). It is shown that the fluctuations of 250- and 340- loops in MT-II system were similar to that of MT systems in 20 ns simulations (Figure S12). Three MD trajectories of H274Y systems (MT, WT→MT, and MT-II)


30

Page 31 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


showed the high fluctuations of 250- and 340- loops, which was not observed in two WT trajectories (WT and HY→WT). Therefore, H274Y mutation might be the source of these distinctive behaviors. The calculated RMSF values were comparable with the experimental B-factor of the crystal structures 3CKZ and 3CL0 (WT) as well as 2HTY and 2HU4 (MT) which also showed that 250loop in H274Y’s crystal structures were higher than WT (Figure S13). As mentioned above, several unusual peptide bonds in 340-loop’s region in H274Y crystal structures21 (Figure S13ab), also implied higher motions in MT than in WT (Figure S13c-d). OTV was designed to improve the drug bioavailability74 since the pioneer drug ZMR was too polar despite its high efficacy. Unfortunately, despite the high affinity of OTV on WT, it became vulnerable to drug resistance due to mutations in NA. The present study also suggests the importance of glycerol-like moiety (or any appropriate H-bond acceptor) at the C6 position in substrate-like compounds for the NA inhibitor, especially to deal with H274Y mutant. This was demonstrated by docking seventeen OTV analogues75 (Table S3) to the Box-1 and Box-2 of the last snapshot of the MT system. Docking results showed that unlike OTV, these derivatives (except OT2, OT3 and OT15) were found to bind at the entrance of Path-1. Moreover, their hydroxyl group(s) at C6 position formed H-bond with R371 (Figure S14). Recently, laninamivir (LNV)76,77, a potent NA inhibitor, was designed with glycerol-like moiety, replacing one of its hydroxyl with a methoxy group. In an enzymatic assay study, LNV was proven to be still effective toward the H274Y mutant78. There was one H-bond between the hydroxyl group of LNV and R371 when LNV docked to the MT snapshot at 20 ns (Figure S15). CONCLUSION


31


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 55

This study aimed to investigate the mechanism of OTV resistance, starting from the simulation of NA in apo form prior to drug binding. Initially, the role of native-H274 in the NA binding site was studied using MD simulations of WT systems in all possible protonation states of histidines. It was shown that the optimum stability of the NA structure was achieved when all histidines were in neutral states with protonated-εN. Specifically, H274 played an important role in maintaining the optimum H-bond network among the 250, 270 and 300 loops. Furthermore, MD simulation of the H274Y system (MT) using the holo form of the NA crystal structure showed several distinctive behaviors of the mutant. These behaviors included the high fluctuation of loops around Y274 and the H-bond formation of E276-R292. The other two MD systems (WT→MT and HY→WT) were simulated and compared to WT and MT to ascertain that these distinct behaviors were due to the H274Y mutation and not due to the removal of the ligand from the holo form of NA. H-bond analysis, dihedral analysis, and essential dynamics studies were performed on the MD trajectories of these four systems (WT, MT, WT→MT and HY→WT) to understand the direct effects of substituting histidine with tyrosine at this position on the integrity of the NA binding site. Sliding box docking was performed to correlate the mutation effects toward OTV binding. The results demonstrated the disruption of the NA binding site in apo form. The disruption was initiated by the repulsive effect of Y274 on the 250-loop, which led to the changes in the H-bond network around residue 274. Noticeable changes on the rotation of the W295 side chain and upward movement of the 340-loop were also observed, phenomena that were not reported previously. Moreover, the sliding box docking results demonstrated that the predicted binding pathway of OTV, but not ZMR, was compromised in the H274Y system. The significance of the functional group at position C6 in OTV resistance was highlighted. This work also underlined the importance of studying the direct effect of mutation on the dynamics of the


32

Page 33 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


protein itself prior to drug binding as a complement to the study of ligand-receptor interaction at the active site. It is hoped that these results could improve the understanding of OTV resistances and shed some light in the design of novel anti-influenza drug. ASSOCIATED CONTENT Supporting Information Figure S1: The grid-box setup in sliding box docking. Figure S2: Scheme of NA binding site, focused at loops and residues around the H274Y point of mutation. Figure S3: The upward movement of the 340-loop in the MT system. Figure S4-S5: Sliding box docking results of OTV in WT and MT systems, respectively. Figure S6: Plot of docking energies of OTV binding pathways in WT. Figure S7: Sliding box docking results of PMV into MT system. Figure S8: Inspection of the 340-loop from coordinates and electron density in the crystal structure. Figure S9: Overall mechanism of the NA binding site disruption due to H274Y. Figure S10: Extended conformation of WT→MT at 40 ns. Figure S11: The conformation of 340-loop in MT system throughout 30 ns of simulation. Figure S12: Relative comparison of B-factors of crystal structures of WT and H274Y. Figure S13: RMSF of a separate 63 ns simulation of MT system. Figure S14: Docking results OTV derivatives at the entrance of NA active site. Figure S15. Predicted binding pathway of LNV to the MT system at 20 ns. Table S1: H-bond analysis. Table S2: Pairwise decomposition of interaction energies. Table S3: List of OTV derivatives that docked to the entrance of Path-1. Movie S1: The timeline mechanism of H274Y effects to the surrounding residues. These materials are available free of charge via the internet at http://pubs.acs.org.


33


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 55

AUTHOR INFORMATION Corresponding Author *Habibah A. Wahab. Email: [email protected]; [email protected]. Funding Sources Malaysian Ministry of Science, Technology and Innovation and Universiti Sains Malaysia. ACKNOWLEDGMENT H.A.W. and N. M. gratefully acknowledges the following grants by Malaysian Ministry of Science, Technology and Innovation and Universiti Sains Malaysia which have been used to disburse GRA/honoraria for M.Y. throughout his non-exchange programme studentship as well as PhD study in USM: (304/PFARMASI/650497/I126 (from 1 April 2010 -31 July 2012; 1001/PFARMASI/815091 from 1 Aug 2012-31 March 2013; 304/PFARMASI/650653/I121 from July -31 Dec 2013; 304/PFARMASI/650545/I121 Mar- June 2014; and others during the period of his study). We would also like to thank Drs Vannajan Lee, Abdul Mutalib, and Ukun MS Soedjanaatmadja for their helpful suggestions. REFERENCES 1. Taubenberger, J.; Morens, D., 1918 Influenza: the Mother of All Pandemics. Emerg. Infect. Dis. 2006, 12, 15-22. 2. HHS H5N1 Avian Flu (H5N1 Bird Flu). http://www.flu.gov/about_the_flu/h5n1/ (July 15, 2013), 3. WHO Background and Summary of Human Infection with Influenza A(H7N9) Virus– as of 5 April 2013. http://www.who.int/influenza/human_animal_interface/update_20130405/en/index.html (June 27, 2013), 4. CDC Avian Influenza A (H7N9) Virus. http://www.cdc.gov/flu/avianflu/h7n9-virus.htm (July 15, 2013), 5. WHO Number of Confirmed Human Cases of Avian Influenza A(H7N9) Reported to WHO. http://www.who.int/influenza/human_animal_interface/influenza_h7n9/Data_Reports/en/index.h tml (June 27, 2013),


34

Page 35 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6. Itzstein, M.; Thomson, R. Anti-Influenza Drugs: The Development of Sialidase Inhibitors. In Antiviral Strategies, Kräusslich, H.-G.; Bartenschlager, R., Eds.; Springer Berlin Heidelberg: 2009; Vol. 189, Chapter 5, pp 111-154. 7. Tong, S.; Zhu, X.; Li, Y.; Shi, M.; Zhang, J.; Bourgeois, M.; Yang, H.; Chen, X.; Recuenco, S.; Gomez, J.; Chen, L.-M.; Johnson, A.; Tao, Y.; Dreyfus, C.; Yu, W.; McBride, R.; Carney, P. J.; Gilbert, A. T.; Chang, J.; Guo, Z.; Davis, C. T.; Paulson, J. C.; Stevens, J.; Rupprecht, C. E.; Holmes, E. C.; Wilson, I. A.; Donis, R. O., New World Bats Harbor Diverse Influenza A Viruses. PLoS Pathog. 2013, 9, e1003657. 8. Matsubara, T.; Onishi, A.; Saito, T.; Shimada, A.; Inoue, H.; Taki, T.; Nagata, K.; Okahata, Y.; Sato, T., Sialic Acid-Mimic Peptides As Hemagglutinin Inhibitors for AntiInfluenza Therapy. J. Med. Chem. 2010, 53, 4441-4449. 9. Gambaryan, A.; Tuzikov, A.; Pazynina, G.; Desheva, J.; Bovin, N.; Matrosovich, M.; Klimov, A., 6-sulfo Sialyl Lewis X Is the Common Receptor Determinant Recognized by H5, H6, H7 and H9 Influenza Viruses of Terrestrial Poultry. Virol. J. 2008, 5, 85. 10. Yusuf, M.; Konc, J.; Sy Bing, C.; Trykowska Konc, J.; Ahmad Khairudin, N. B.; Janezic, D.; Wahab, H. A., Structurally Conserved Binding Sites of Hemagglutinin as Targets for Influenza Drug and Vaccine Development. J. Chem. Inf. Model. 2013, 53, 2423-2436. 11. Hayden, F. G., Antiviral Resistance in Influenza Viruses — Implications for Management and Pandemic Response. N. Engl. J. Med. 2006, 354, 785-788. 12. Moscona, A., Oseltamivir Resistance — Disabling Our Influenza Defenses. N. Engl. J. Med. 2005, 353, 2633-2636. 13. Soundararajan, V.; Tharakaraman, K.; Raman, R.; Raguram, S.; Shriver, Z.; Sasisekharan, V.; Sasisekharan, R., Extrapolating from Sequence—the 2009 H1N1 'swine' Influenza Virus. Nat. Biotech. 2009, 27, 510-513. 14. Hurt, A. C.; Ernest, J.; Deng, Y.-M.; Iannello, P.; Besselaar, T. G.; Birch, C.; Buchy, P.; Chittaganpitch, M.; Chiu, S.-C.; Dwyer, D.; Guigon, A.; Harrower, B.; Kei, I. P.; Kok, T.; Lin, C.; McPhie, K.; Mohd, A.; Olveda, R.; Panayotou, T.; Rawlinson, W.; Scott, L.; Smith, D.; D'Souza, H.; Komadina, N.; Shaw, R.; Kelso, A.; Barr, I. G., Emergence and Spread of Oseltamivir-resistant a(H1N1) Influenza Viruses in Oceania, South East Asia and South Africa. Antiviral Res. 2009, 83, 90-93. 15. Liu, X.; Li, T.; Zheng, Y.; Wong, K.-W.; Lu, S.; Lu, H., Poor Responses to Oseltamivir Treatment in a Patient with Influenza a (H7N9) Virus Infection. Emerg Microbes Infect 2013, 2, e27. 16. Abed, Y.; Baz, M.; Boivin, G., Impact of Neuraminidase Mutations Conferring Influenza Resistance to Neuraminidase Inhibitors in the N1 and N2 Genetic Backgrounds. Antivir. Ther. 2006, 11, 971-976. 17. Shan, Y.; Ma, Y.; Wang, M.; Dong, Y., Recent Advances in the Structure-Based Design of Neuraminidase Inhibitors as Antiinfluenza Agents. Curr Med Chem 2012, 19, 5885-5894. 18. Ikram, N. K. K.; Durrant, J. D.; Muchtaridi, M.; Zalaludin, A. S.; Purwitasari, N.; Mohamed, N.; Rahim, A. S. A.; Lam, C. K.; Normi, Y. M.; Rahman, N. A.; Amaro, R. E.; Wahab, H. A., A Virtual Screening Approach For Identifying Plants with Anti H5N1 Neuraminidase Activity. J. Chem. Inf. Model. 2015, 55, 308-316. 19. Frimayanti, N.; Lee, V.; Zain, S.; Wahab, H.; Abd. Rahman, N., 2D, 3D-QSAR, and Pharmacophore Studies on Thiazolidine-4-carboxylic Acid Derivatives as Neuraminidase Inhibitors in H3N2 Influenza Virus. Med. Chem. Res. 2014, 23, 1447-1453.


35


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 55

20. Wang, M. Z.; Tai, C. Y.; Mendel, D. B., Mechanism by Which Mutations at His274 Alter Sensitivity of Influenza A Virus N1 Neuraminidase to Oseltamivir Carboxylate and Zanamivir. Antimicrob. Agents Chemother. 2002, 46, 3809-3816. 21. Collins, P. J.; Haire, L. F.; Lin, Y. P.; Liu, J.; Russell, R. J.; Walker, P. A.; Skehel, J. J.; Martin, S. R.; Hay, A. J.; Gamblin, S. J., Crystal Structures of Oseltamivir-resistant Influenza Virus Neuraminidase Mutants. Nature 2008, 453, 1258-1261. 22. Rungrotmongkol, T.; Malaisree, M.; Nunthaboot, N.; Sompornpisut, P.; Hannongbua, S., Molecular Prediction of Oseltamivir Efficiency Against Probable Influenza a (H1N1-2009) Mutants: Molecular Modeling Approach. Amino Acids 2010, 39, 393-398. 23. Malaisree, M.; Rungrotmongkol, T.; Nunthaboot, N.; Aruksakunwong, O.; Intharathep, P.; Decha, P.; Sompornpisut, P.; Hannongbua, S., Source of Oseltamivir Resistance in Avian Influenza H5N1 Virus with the H274Y Mutation. Amino Acids 2009, 37, 725-732. 24. Le, L.; Lee, E. H.; Hardy, D. J.; Truong, T. N.; Schulten, K., Molecular Dynamics Simulations Suggest that Electrostatic Funnel Directs Binding of Tamiflu to Influenza N1 Neuraminidases. PLoS Comput. Biol. 2010, 6, e1000939. 25. Karthick, V.; Shanthi, V.; Rajasekaran, R.; Ramanathan, K., Exploring the Cause of Oseltamivir Resistance Against Mutant H274Y Neuraminidase by Molecular Simulation Approach. Appl Biochem Biotechnol 2012, 167, 237-249. 26. Wang, N. X.; Zheng, J. J., Computational Studies of H5N1 Influenza Virus Resistance to Oseltamivir. Protein Sci. 2009, 18, 707-715. 27. Park, J. W.; Jo, W. H., Infiltration of Water Molecules into the Oseltamivir-Binding Site of H274Y Neuraminidase Mutant Causes Resistance to Oseltamivir. J. Chem. Inf. Model. 2009, 49, 2735-2741. 28. Kati, W. M.; Saldivar, A. S.; Mohamadi, F.; Sham, H. L.; Laver, W. G.; Kohlbrenner, W. E., GS4071 Is a Slow-Binding Inhibitor of Influenza Neuraminidase from Both A and B Strains. Biochem. Biophys. Res. Commun. 1998, 244, 408-413. 29. Li, L.; Li, Y.; Zhang, L.; Hou, T., Theoretical Studies on the Susceptibility of Oseltamivir against Variants of 2009 A/H1N1 Influenza Neuraminidase. J. Chem. Inf. Model. 2012, 52, 2715-2729. 30. Ripoll, D. R.; Khavrutskii, I. V.; Chaudhury, S.; Liu, J.; Kuschner, R. A.; Wallqvist, A.; Reifman, J., Quantitative Predictions of Binding Free Energy Changes in Drug-Resistant Influenza Neuraminidase. PLoS Comput. Biol. 2012, 8, e1002665. 31. Amaro, R. E.; Cheng, X.; Ivanov, I.; Xu, D.; McCammon, J. A., Characterizing Loop Dynamics and Ligand Recognition in Human- and Avian-Type Influenza Neuraminidases via Generalized Born Molecular Dynamics and End-Point Free Energy Calculations. J. Am. Chem. Soc. 2009, 131, 4702-4709. 32. Sung, J. C.; Wynsberghe, A. W. V.; Amaro, R. E.; Li, W. W.; McCammon, J. A., Role of Secondary Sialic Acid Binding Sites in Influenza N1 Neuraminidase. J. Am. Chem. Soc. 2010, 132, 2883-2885. 33. Cheng, L. S.; Amaro, R. E.; Xu, D.; Li, W. W.; Arzberger, P. W.; McCammon, J. A., Ensemble-Based Virtual Screening Reveals Potential Novel Antiviral Compounds for Avian Influenza Neuraminidase. J. Med. Chem. 2008, 51, 3878-3894. 34. Gubareva, L. V.; Robinson, M. J.; Bethell, R. C.; Webster, R. G., Catalytic and Framework Mutations in the Neuraminidase Active Site of Influenza Viruses That Are Resistant to 4-guanidino-neu5ac2en. J. Virol. 1997, 71, 3385-3390.


36

Page 37 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


35. Klenk, H. D.; Matrosovich, M. N.; Stech, J., Avian Influenza. Karger: Basel, Switzerland, 2008. 36. Woods, C. J.; Malaisree, M.; Long, B.; McIntosh-Smith, S.; Mulholland, A. J., Analysis and Assay of Oseltamivir-Resistant Mutants of Influenza Neuraminidase via Direct Observation of Drug Unbinding and Rebinding in Simulation. Biochemistry 2013, 52, 8150-8164. 37. Woods, C. J.; Malaisree, M.; Pattarapongdilok, N.; Sompornpisut, P.; Hannongbua, S.; Mulholland, A. J., Long Time Scale GPU Dynamics Reveal the Mechanism of Drug Resistance of the Dual Mutant I223R/H275Y Neuraminidase from H1N1-2009 Influenza Virus. Biochemistry 2012, 51, 4364-4375. 38. Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E., The Protein Data Bank. Nucleic Acids Res. 2000, 28, 235-242. 39. Russell, R. J.; Haire, L. F.; Stevens, D. J.; Collins, P. J.; Lin, Y. P.; Blackburn, G. M.; Hay, A. J.; Gamblin, S. J.; Skehel, J. J., The Structure of H5N1 Avian Influenza Neuraminidase Suggests New Opportunities for Drug Design. Nature 2006, 443, 45-49. 40. Accelrys Discovery Studio, 2.5; Accelrys Software Inc.: San Diego, USA, 2007. 41. 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.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Kolossváry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; J., L.; 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.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER, 11; University of California, San Francisco, 2010. 42. Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P., A Point-charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-phase Quantum Mechanical Calculations. J Comput Chem. 2003, 24, 1999-2012. 43. Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C., Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-alkanes. J. Comput. Phys. 1977, 23, 327-341. 44. Loncharich, R. J.; Brooks, B. R.; Pastor, R. W., Langevin Dynamics of Peptides: the Frictional Dependence of Isomerization Rates of N-acetylalanyl-n2-methylamide. Biopolymers 1992, 32, 523-535. 45. Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R., Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684-3690. 46. Darden, T.; York, D.; Pedersen, L., Particle Mesh Ewald: An N⊕log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089-10092. 47. Miller, B. R.; McGee, T. D.; Swails, J. M.; Homeyer, N.; Gohlke, H.; Roitberg, A. E., MMPBSA.py: An Efficient Program for End-State Free Energy Calculations. J. Chem. Theory Comput 2012, 8, 3314-3321. 48. Bakan, A.; Meireles, L. M.; Bahar, I., ProDy: Protein Dynamics Inferred from Theory and Experiments. Bioinformatics 2011, 27, 1575-1577. 49. Humphrey, W.; Dalke, A.; Schulten, K., VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33-38. 50. Tran, D.-T.; Le, L.; Truong, T., Discover Binding Pathways Using the Sliding Bindingbox Docking Approach: Application to Binding Pathways of Oseltamivir to Avian Influenza H5N1 Neuraminidase. J. Comput. Aided Mol. Des. 2013, 27, 689-695.


37


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 55

51. Trott, O.; Olson, A. J., Autodock Vina: Improving the Speed and Accuracy of Docking with a New Scoring Function, Efficient Optimization, and Multithreading. J. Comp. Chem. 2010, 31, 455-461. 52. Accelrys Discovery Studio Visualizer, 4.1; Accelrys Software Inc.: San Diego, USA, 2014. 53. Uranga, J.; Mikulskis, P.; Genheden, S.; Ryde, U., Can the Protonation State of Histidine Residues Be Determined from Molecular Dynamics Simulations? Computational and Theoretical Chemistry 2012, 1000, 75-84. 54. Amaro, R. E.; Minh, D. D. L.; Cheng, L. S.; Lindstrom, W. M.; Olson, A. J.; Lin, J.-H.; Li, W. W.; McCammon, J. A., Remarkable Loop Flexibility in Avian Influenza N1 and Its Implications for Antiviral Drug Design. J. Am. Chem. Soc. 2007, 129, 7764-7765. 55. Han, N.; Mu, Y., Plasticity of 150-Loop in Influenza Neuraminidase Explored by Hamiltonian Replica Exchange Molecular Dynamics Simulations. PLoS One 2013, 8, e60995. 56. Dolinsky, T. J.; Nielsen, J. E.; McCammon, J. A.; Baker, N. A., PDB2PQR: an Automated Pipeline for the Setup of Poisson–boltzmann Electrostatics Calculations. Nucleic Acids Res. 2004, 32, W665-W667. 57. Søndergaard, C. R.; Olsson, M. H. M.; Rostkowski, M.; Jensen, J. H., Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values. J. Chem. Theory Comput 2011, 7, 2284-2295. 58. Olsson, M. H. M.; Søndergaard, C. R.; Rostkowski, M.; Jensen, J. H., PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions. J. Chem. Theory Comput 2011, 7, 525-537. 59. Amaro, R.; Li, W. Molecular-Level Simulation of Pandemic Influenza Glycoproteins. In Computational Drug Discovery and Design, Baron, R., Ed.; Springer New York: 2012; Vol. 819, Chapter 34, pp 575-594. 60. Wu, Y.; Qin, G.; Gao, F.; Liu, Y.; Vavricka, C. J.; Qi, J.; Jiang, H.; Yu, K.; Gao, G. F., Induced Opening of Influenza Virus Neuraminidase N2 150-loop Suggests an Important Role in Inhibitor Binding. Sci. Rep. 2013, 3. 61. Li, Q.; Qi, J.; Wu, Y.; Kiyota, H.; Tanaka, K.; Suhara, Y.; Ohrui, H.; Suzuki, Y.; Vavricka, C. J.; Gao, G. F., Functional and Structural Analysis of Influenza Neuraminidase N3 Offers Further Insight into the Mechanisms of Oseltamivir-resistance. J. Virol. 2013, 87, 1001610024. 62. Bissantz, C.; Kuhn, B.; Stahl, M., A Medicinal Chemist’s Guide to Molecular Interactions. J. Med. Chem. 2010, 53, 5061-5084. 63. Mitchell, J. B. O.; Nandi, C. L.; McDonald, I. K.; Thornton, J. M.; Price, S. L., Amino/Aromatic Interactions in Proteins: Is the Evidence Stacked Against Hydrogen Bonding? J. Mol. Biol. 1994, 239, 315-331. 64. Udommaneethanakit, T.; Rungrotmongkol, T.; Bren, U.; Frecer, V.; Stanislav, M., Dynamic Behavior of Avian Influenza A Virus Neuraminidase Subtype H5N1 in Complex with Oseltamivir, Zanamivir, Peramivir, and Their Phosphonate Analogues. J. Chem. Inf. Model. 2009, 49, 2323-2332. 65. Stoll, V.; Stewart, K. D.; Maring, C. J.; Muchmore, S.; Giranda, V.; Gu, Y.-g. Y.; Wang, G.; Chen, Y.; Sun, M.; Zhao, C.; Kennedy, A. L.; Madigan, D. L.; Xu, Y.; Saldivar, A.; Kati, W.; Laver, G.; Sowin, T.; Sham, H. L.; Greer, J.; Kempf, D., Influenza Neuraminidase Inhibitors: Structure-Based Design of a Novel Inhibitor Series. Biochemistry 2003, 42, 718-727.


38

Page 39 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


66. Martins-José, A., Sliding Box Docking: a New Stand-alone Tool for Managing Dockingbased Virtual Screening along the DNA Helix Axis. Bioinformation 2013, 9, 750-751. 67. Abed, Y.; Baz, M.; Boivin, G., A Novel Neuraminidase Deletion Mutation Conferring Resistance to Oseltamivir in Clinical Influenza A/H3N2 Virus. J. Infect. Dis. 2009, 199, 180183. 68. Hay, A.; Collins, P.; Russell, R. Antivirals and Resistance. In Avian Influenza, Klenk, H. D., Ed.; Karger: Basel, Switzerland, 2008; Vol. 27, pp 252-271 69. Smith, B. J.; Huyton, T.; Joosten, R. P.; McKimm-Breschkin, J. L.; Zhang, J.-G.; Luo, C. S.; Lou, M.-Z.; Labrou, N. E.; Garrett, T. P. J., Structure of a Calcium-deficient Form of Influenza Virus Neuraminidase: Implications for Substrate Binding. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2006, 62, 947-952. 70. Lawrenz, M.; Wereszczynski, J.; Amaro, R.; Walker, R.; Roitberg, A.; McCammon, J. A., Impact of Calcium on N1 Influenza Neuraminidase Dynamics and Binding Free Energy. Proteins: Struct., Funct., Bioinf. 2010, 78, 2523-2532. 71. Berg, J. M.; Tymoczko, J. L.; Stryer, L. Section 3.2, Primary Structure: Amino Acids Are Linked by Peptide Bonds to Form Polypeptide Chains. In Biochemistry; W H Freeman: New York, 2002. 72. Graf, M. H.; Bren, U.; Haltrich, D.; Oostenbrink, C., Molecular Dynamics Simulations Give Insight into D-glucose Dioxidation at C2 and C3 by Agaricus Meleagris Pyranose Dehydrogenase. J. Comput. Aided Mol. Des. 2013, 27, 295-304. 73. Bren, U.; Oostenbrink, C., Cytochrome P450 3A4 Inhibition by Ketoconazole: Tackling the Problem of Ligand Cooperativity Using Molecular Dynamics Simulations and Free-Energy Calculations. J. Chem. Inf. Model. 2012, 52, 1573-1582. 74. Kim, C. U.; Lew, W.; Williams, M. A.; Liu, H.; Zhang, L.; Swaminathan, S.; Bischofberger, N.; Chen, M. S.; Mendel, D. B.; Tai, C. Y.; Laver, W. G.; Stevens, R. C., Influenza Neuraminidase Inhibitors Possessing a Novel Hydrophobic Interaction in the Enzyme Active Site: Design, Synthesis, and Structural Analysis of Carbocyclic Sialic Acid Analogues with Potent Anti-Influenza Activity. J. Am. Chem. Soc. 1997, 119, 681-690. 75. Du, Q.-S.; Wang, S.-Q.; Chou, K.-C., Analogue Inhibitors by Modifying Oseltamivir Based on the Crystal Neuraminidase Structure for Treating Drug-resistant H5N1 Virus. Biochem. Biophys. Res. Commun. 2007, 362, 525-531. 76. Yamashita, M.; Tomozawa, T.; Kakuta, M.; Tokumitsu, A.; Nasu, H.; Kubo, S., CS8958, a Prodrug of the New Neuraminidase Inhibitor R-125489, Shows Long-Acting AntiInfluenza Virus Activity. Antimicrob. Agents Chemother. 2009, 53, 186-192. 77. Yamashita, M., Laninamivir and Its Prodrug, CS-8958: Long-acting Neuraminidase Inhibitors for the Treatment of Influenza. Antiviral Chem. Chemother. 2010, 21, 71-84. 78. Vavricka, C. J.; Li, Q.; Wu, Y.; Qi, J.; Wang, M.; Liu, Y.; Gao, F.; Liu, J.; Feng, E.; He, J.; Wang, J.; Liu, H.; Jiang, H.; Gao, G. F., Structural and Functional Analysis of Laninamivir and its Octanoate Prodrug Reveals Group Specific Mechanisms for Influenza NA Inhibition. PLoS Pathog. 2011, 7, e1002249.


39


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 55

TABLES AND FIGURES Table 1. List of MD systems to elucidate the mechanism of H274Y in apo state of NA.

1.

HID

WT, all histidines protonated at δN

PDB ID 2HTY

2.

HIE or WT

WT, all histidines protonated at εN

2HTY

20

3.

HIP

WT, all histidines protonated at δN and εN

2HTY

20

4.

MT

Crystal structure of H274Y in holo form

3CKZ

20

5.

WT→MT

WT mutated into H274Y

2HTY

20

6.

HY→WT

MT mutated back into WT (Y274H)

3CKZ

20

No.

System

Explanation

Timescale (ns)* 20

* 20 ns of production run, in addition to 150 ps of heating and 1 ns of equilibration stages

Table 2. Time evolution of hydrogen bond formation among residues within 5 Å from position 274 in WT, MT, WT→MT, and HY→WT systems throughout 20 ns of MD simulation, which shows significant differences between the WT and H274Y systems.


40

Page 41 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


System Acceptor** WT HY→WT MT WT→MT WT HY→WT MT WT→MT WT HY→WT MT WT→MT WT HY→WT MT WT→MT WT HY→WT MT WT→MT WT HY→WT MT WT→MT WT HY→WT MT WT→MT WT HY→WT MT WT→MT

T242@OG1 T242@OG1 T242@OG1 T242@OG1 S246@OG S246@OG S246@OG S246@OG S246@OG S246@OG S246@OG S246@OG Y252@OH Y252@OH Y252@OH Y252@OH R224@NH1 R224@NH1 R224@NH1 R224@NH1 R224@NH1 R224@NH1 R224@NH1 R224@NH1 R224@NH1 R224@NH1 R224@NH1 R224@NH1 R224@NH1 R224@NH1 N247@ND2 N247@ND2 N247@ND2 N247@ND2 N272@OD1 N272@OD1 N272@OD1 N272@OD1

Donor** Y275@O Y275@O Y275@O Y275@O E276@OE2 E276@OE1 E276@OE2 E276@OE1 E276@OE2 E276@OE1 E276@OE1 E276@OE2 Y273@O Y273@O Y273@O Y273@O H274@NE2 H274@NE2 Y274@OH Y274@OH P245@O P245@O P245@O P245@O S246@OG S246@OG S246@OG S246@OG S246@OG S246@OG N294@O N294@O N294@O N294@O Y274@OH Y274@OH Y274@OH Y274@OH

% Occ* 64.8 65.0 18.6 51.4 58.8 40.3 38.6 30.7 9.5 4.2 16.8 10.8 10.9 30.0 3.8 4.2 0.0 0.0 2.9 9.7 10.7 17.2 0.0 1.5 1.9 1.1 4.6 2.1 3.9 2.8 25.8 10.2 0.0 0.0 0.0 0.0 11.7 0.0

2

4

6

8

Timescale (ns) 10 12

14

16

18

20

Scale 1234567 * % Occ = percentage of occupancy; Scale 1 = 0-5%; 2 = 5-20%; 3 = 20-40%; 4 = 40-60%; 5 = 60-80%; 6 = 80-95%; 7 = 95-100% occupancy ** Members of the 250, 270, and 300 loops are colored in orange, red, and purple respectively. Other residues located around the 250-loop are colored in yellow.

Table 2. Continued.


41


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

System Acceptor**

Donor**

H274@NE2 E276@OE2 WT HY→WT H274@NE2 E276@OE2 H274@NE2 E276@OE1 Y274@OH E276@OE2 MT WT→MT Y274@OH E276@OE1 Y274@OH E276@OE2 N272@O H296@NE2 WT HY→WT N272@O H296@NE2 N272@OD1 H296@NE2 N272@O H296@NE2 MT WT→MT N272@O H296@NE2 E276@OE2 N294@ND2 WT E276@OE1 N294@ND2 HY→WT E276@OE1 N294@ND2 E276@OE2 N294@ND2 E276@OE2 N294@ND2 MT E276@OE1 N294@ND2 WT→MT E276@OE2 N294@ND2 E276@OE1 N294@ND2 E276@OE1 N294@ND2 E276@OE2 N294@ND2 E276@OE2 R292@NH2 WT HY→WT E276@OE2 R292@NH2 E276@OE2 R292@NH2 MT E276@OE2 R292@NH1 E276@OE1 R292@NH1 WT→MT E276@OE2 R292@NH2 E276@OE2 R292@NH1 E276@OE1 R292@NH1 E276@OE1 R292@NH2 R292@NE N294@OD1 WT R292@NH1 N294@OD1 HY→WT R292@NH1 N294@OD1 R292@NE N294@OD1 R292@NH1 N294@OD1 MT R292@NH2 N294@OD1 R292@NH2 N294@ND2 WT→MT R292@NH1 N294@OD1 R292@NE N294@OD1 R292@NE D293@O WT HY→WT R292@NE D293@O R292@NH2 D293@O MT R292@NE D293@O WT→MT R292@NE D293@O R292@NH1 D293@O

% Occ* 0.0 41.9 24.4 0.0 49.5 46.0 44.5 58.7 3.8 0.0 0.0 4.7 2.7 12.4 1.9 47.3 9.0 6.9 6.4 3.1 2.1 0.0 0.0 32.2 4.1 2.2 6.8 6.3 5.7 4.6 35.7 30.6 57.7 51.9 5.3 4.0 1.0 18.1 10.2 0.0 0.0 58.1 7.3 17.4 3.0

2

4

6

8

Timescale (ns) 10 12

Page 42 of 55

14

16

18

20

Scale 1234567 * % Occ = percentage of occupancy; Scale 1 = 0-5%; 2 = 5-20%; 3 = 20-40%; 4 = 40-60%; 5 = 60-80%; 6 = 80-95%; 7 = 95-100% occupancy ** Members of the 250, 270, and 300 loops are colored in orange, red, and purple respectively. Other residues located around the 250-loop are colored in yellow.


42

Page 43 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Binding mode of interaction of OTV (green colored sticks) on the active site of NA and the overview of the hypotheses of H274Y effect towards the mechanism of OTV resistance.

Figure 2. Four possible interactions of H274 with neighboring residues from the 250-loop and the 300-loop (visualized in green colored lines) based on the crystal structure of PDB ID 2HU439.


43


2

RMSD (Å)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 55

1.5 1 0.5 0 0

2

4

6

8

10

12

14

16

18

20

Time (ns) 100 per. Mov. Avg. (HID)

100 per. Mov. Avg. (HIE)

100 per. Mov. Avg. (HIP)

100 per. Mov. Avg. (MT)

Figure 3. Time evolution plot of RMSD protein backbone throughout 20 ns of MD simulation. The bold lines are the moving averages over 100 data points.

Figure 4. RMSF of backbone for each residue of NA throughout 20 ns of simulation. The marked differences in RMSF are highlighted using boxes.


44

Page 45 of 55

Figure 5. The average structures of (a) HID, (b) HIE, (c) HIP, and (d) MT systems and the time evolution snapshots of (e) HID, (f) HIE, (g) HIP, and (h) MT throughout 20 ns of simulations. Hydrogen bonds and salt bridges are depicted in green dashes and yellow lines, respectively.

2 1.5

RMSD (Å)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 0.5 0 0

2

4

6

8

100 per. Mov. Avg. (WT) 100 per. Mov. Avg. (MT)

10 12 14 16 18 20 Time (ns) 100 per. Mov. Avg. (WT→MT) 100 per. Mov. Avg. (HY→WT)

Figure 6. The time evolution RMSD of NA backbone in four studied systems: WT, MT, WT→MT, and HY→WT. The bold lines are moving averages over 100 data points.


45


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 55

Figure 7. RMSF values of protein backbone in WT (green color), WT→MT (red color), MT (orange color), and HY→WT (blue color) systems. Loops around position 274, (i.e., the 250, 270, 300, and 340 loops) are highlighted in separate boxes.


46

Page 47 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8. The superimposed conformations of (a) WT, (b) WT→MT, (c) MT, and (d) HY→WT extracted from four ns intervals. The time evolution snapshots are represented by gradient redwhite-blue colors.


47


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 48 of 55

Figure 9. Illustration of OTV binding with (a) WT and (b) final snapshot of MT system. The Hbond of E276-R292 formed polar environment at hydrophobic subpocket for OTV.


48

Page 49 of 55

WT

Energy (kcal/mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


HY→WT

MT

WT→MT

1.0 0.0 -1.0 -2.0 -3.0 -4.0 -5.0 -6.0 -7.0 -8.0 -9.0 276*-292

224*-245

295*-296

294**-247

292*-293

242*-275

294*-292

* side-chain; ** backbone

Figure 10. Pairwise decomposition of the interaction energies between two residues.


49


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 50 of 55

Figure 11. Dihedral angle plot of (a) H296, (b) W295, and (c) H/Y274 over 20 ns of simulation. The time evolution conformational changes of these residues and loops in (d) WT→MT and (e) MT systems are presented using gradient colors red-white-blue.


50

Page 51 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 12. Essential dynamics of (a) WT, (b) MT, (c) WT→MT, and (d) HY→WT systems. Area inside dashed box covers the loops around position 274. The arrows represent the amplitude and direction of movement.


51


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 52 of 55

Figure 13. Plot of sliding-docking results: (a) OTV and (b) ZMR on WT system at 0 ns; (c) OTV and (d) ZMR on WT system at 20 ns; (e) OTV and (f) ZMR on MT system at 20 ns. Boxes present the empty binding site surface, which are composed of the 150, 250, 270, 300, and 340 loops. Binding pathways are represented by the gradient color of ligands: red to orange to green colors.


52

Page 53 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 14. The role of ZMR glycerol moiety in the interaction with R371 to reach the active site through the 340-loop in MT system.


53


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 54 of 55

(a)

(b)

Figure 15. Endpoint interaction of (a) OTV and (b) ZMR with the last snapshot of the MT system at 20 ns. The hydrogen bond is represented by green dashed lines. Docked poses of OTV and ZMR are visualized in green colored sticks, and the co-crystal structure of OTV is visualized in gray.


54

Page 55 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For Table of Contents Use Only

H274Y’s effect on Oseltamivir Resistance: What Happens Before the Drug Enters the Binding Site Muhammad Yusuf, Nornisah Mohamed, Suriyati Mohamad, Dusanka Janezic, K.V. Damodaran and Habibah A. Wahab*

The direct effect of the H274Y mutation to the dynamics behavior of the NA binding site was investigated using molecular dynamics simulation. Several changes in the mutant system were observed as compared to the wild type, including changes in the H-bond networks and loops motion. These structural changes were suggested to disrupt the binding pathway of OTV.


55

H274Y's Effect on Oseltamivir Resistance: What Happens Before the

Recommend Documents