Understanding the cross-resistance of oseltamivir to H1N1 and H5N1 influenza A neuraminidase mutations using multidimensional computational analyses.

This study embarks on a comprehensive description of the conformational contributions to resistance of neuraminidase (N1) in H1N1 and H5N1 to oseltamivir, using comparative multiple molecular dynamic simulations. The available data with regard to elucidation of the mechanism of resistance as a result of mutations in H1N1 and H5N1 neuraminidases is not well established. Enhanced post-dynamic analysis, such as principal component analysis, solvent accessible surface area, free binding energy calculations, and radius of gyration were performed to gain a precise insight into the binding mode and origin of resistance of oseltamivir in H1N1 and H5N1 mutants. Three significant features reflecting resistance in the presence of mutations H274Y and I222K, of the protein complexed with the inhibitor are: reduced flexibility of the α-carbon backbone; an improved ΔEele of ~15 (kcal/mol) for H1N1 coupled with an increase in ΔGsol (~13 kcal/mol) from wild-type to mutation; a low binding affinity in comparison with the wild-type of ~2 (kcal/mol) and ~7 (kcal/mol) with respect to each mutation for the H5N1 systems; and reduced hydrophobicity of the overall surface structure due to an impaired hydrogen bonding network. We believe the results of this study will ultimately provide a useful insight into the structural landscape of neuraminidase-associated binding of oseltamivir. Furthermore, the results can be used in the design and development of potent inhibitors of neuraminidases.

Introduction

The rapid evolution of highly pathogenic and resistant variants of influenza A viruses, H5N1 and H1N1, takes us to the precipice of a possible new-age influenza pandemic. The first human influenza pandemic documented in the 20th century was the 1918–1919 “Spanish flu”, which claimed the lives of approximately 50 million people worldwide.1–4 Since 1997, the highly pathogenic avian influenza A (H5N1) has had a detrimental effect on the socioeconomic development in countries across Africa, Asia, and the Middle East. Pharmacological advances in viable influenza chemotherapeutics and vaccines could not have predicted the influenza A (H1N1) virus outbreak in 2009–2010, which claimed an estimated 284,500 lives. The viruses H5N1 and H1N1 display similar pathogenesis of respiratory tract infection in humans. However, statistically H5N1 has a higher mortality rate.5

Influenza virus A is highly communicable. H1N1 is easily transmissible between humans, while avian influenza has not yet adapted sufficiently to facilitate human-human spread; the spread of H5N1 is largely dependent on close avian-human contact. Based on the research, the H1N1 and H5N1 viruses can both trace their lineage back to pigs as potential intermediate hosts. Both viruses can undergo reassortment in the intermediate host prior to infection in the new host.6,7 This raises the concern that more virulent and highly transmissible influenza virus strains may infect humans.8

Both the H1N1 and H5N1 viruses express two main surface antigens, H5/H1 (hemagglutinin type 5/1, respectively) and N1 (neuraminidase type 1).9–11 Neuraminidase has been the main target for potential pharmacological interventions.12 It plays a vital role in the spread of infection through cleavage of the viral receptor, sialic acid, from both viral and host proteins. This results in the release of newly replicated virus from an infected cell.12,13 The active site residues of the neuraminidase enzyme are highly conserved within all virus types and subtypes, ie, designed neuraminidase inhibitors may be effective against types A and B of the influenza virus.14–17

There are four effective neuraminidase inhibitors currently available: zanamivir,18 laninamivir,19 peramivir,20 and oseltamivir.21 Zanamivir is commercially known as Relenza® and is administered as an aerosol into the nasal cavity. Zanamivir has been approved in many countries since 1999/2000. Laninamivir is administered in the form of an inhaler and was made available to the Japanese market in 2010 (under the trade name Inavir®). Peramivir, commercially known as PeramiFlu®, is administered by intravenous drip and is approved in Japan, the Republic of Korea, and the People’s Republic of China. Oseltamivir commercially known as Tamiflu® (Figure 1), has long been the “drug of choice” due to ease of oral administration, improved bioavailability, and easy accessibility. There has been an increase in the number of reports describing oseltamivir resistance in H1N1 and H5N1, and the rapid development of resistance toward the drug poses a threat in the event of an outbreak.22,23

Figure 1

Structure of oseltamivir.

Research directed toward gaining an insight into the mechanism of resistance of influenza viruses is imperative in the future design of drugs as well as in the development of pre-emptive measures against possible mutations. The comparative multiple molecular dynamic (MD) simulations conducted in this study offer an atomistic perspective into the complex nature of protein motions as a function of time.24,25 Due to H5N1 virus sharing a 91.47% sequence identity and a conserved binding site with the 2009 H1N1, the results of molecular simulation may apply to both viruses.26

Several computational approaches have been used in an attempt to understand the impact of mutations on oseltamivir resistance to neuraminidase. Karthick et al initiated one of the first MD studies, which was based solely on the natural mutation H274Y of H1N1 isolated during the 2007–2008 Euroasia influenza season.27,28 Malaisree et al refined the source of oseltamivir resistance in avian influenza H5N1 virus with the H274Y mutation using MD simulation.29,30 Wang and Zheng performed a similar study that incorporated the mutations N294S and E119G.31 In the above studies, it was found that although residue 274 is a framework residue in the active site,21 the mutation of histidine (His) to tyrosine (Tyr) significantly impaired the binding of oseltamivir because the binding pocket could not accommodate the steric bulk of the drug. Further to this, the presence of additional mutations was found to amplify resistance of the neuraminidase enzyme to oseltamivir. A large database of crystallized glycoprotein of influenza virus A and B exists, each with a different mutation in the neuraminidase amino acid sequence.26,27,32–34 The MD studies for each mutation revealed enhanced surface expression of oseltamivir-resistant influenza neuraminidase.35 In many instances, these included zanamivir and peramivir resistance.36

Nguyen et al37 and Huang et al38 identified a mutation at position 222 from isoleucine to lysine (Lys) in strains of wild-type (WT) and H274Y mutants of H1N1 and H5N1. The point mutation from isoleucine to Lys exhibited increased resistance toward oseltamivir. Structural residue I222 of H1N1 and H5N1 neuraminidase offered stability to oseltamivir binding.39,40 Structure stability was reinforced by the presence of the hydrophobic pocket created by E276 and R224, which in turn supported a seven hydrogen bond interaction between the active site residues and the drug.30,36

Insufficient information about the origin of resistance of oseltamivir against neuraminidase H1N1 and H5N1 mutations, I222K and H274Y (Figure 2) prompted us to perform a comprehensive multidimensional analysis using a multiple MD approach. Multiple MD simulations have demonstrated better sampling efficiency.41 Multiple trajectories with different initial conditions improved the conformational sampling of MD simulations in proteins.42,43 Several post-dynamic analyses, such as root mean square fluctuation (RMSF), root mean square deviation (RMSD), free binding energy calculations, radius of gyration, solvent accessible surface area (SASA), and principal component analysis (PCA) were performed in order to gain a comprehensive understanding of the impact of mutations on binding and the conformational landscape of the oseltamivir–protein complex. Our research should contribute to understanding the effect of resistance and provide insight into the future development of innovative chemotherapeutics.

Figure 2

Three-dimensional structures of H5N1 and H1N1 neuraminidase A and B, respectively, showing the positions of the studied mutations.

Computational methods

System preparation

The X-ray crystal structures of H1N1 and H5N1 WT and mutant neuraminidase were extracted from Protein Data Bank (PDB) codes 4B7R (WT-H1N1),44 2HU4 (WT-H5N1),15 and 3CL0 (H274Y mutation-H5N1)45 complexed with oseltamivir, obtained from the Research Collaboratory for Structural Bioinformatics. Neuraminidase is a tetrameric enzyme, using the Chimera software package46 a single subunit was selected with the inclusion of an active site drug complex to reduce computational cost. Influenza virus, H1N1 PDB code 4B7R, formed the WT sequence with a His274 residue. Using Chimera, point mutations were introduced at positions 222 (from I to K) and 274 (from H to Y). Table 1 lists the three enzymes with relevant mutations resulting in eight simulations.

Table 1

Crystal structures of the simulated systems, PDB codes, and abbreviations

Simulated system	PDB code	Abbreviation*
H1N1 wild-type	4B7R	WTH1N1
H1N1 I222K	4B7R	I222KH1N1
H1N1 H274Y	4B7R	H274YH1N1
H1N1 H274Y-I222K	4B7R	H274Y-I222KH1N1
H5N1 wild-type	2HU4	WTH5N1
H5N1 I222K	2HU4	I222KH5N1
H5N1 H274Y	3CL0	H274YH5N1
H5N1 H274Y-I222K	3CL0	H274Y-I222KH5N1

Note:

These abbreviations are used throughout the text.

Abbreviation: PDB, Protein Data Bank.

MD simulation

Multiple MD simulations were performed to establish the impact of mutation I222K on the binding of oseltamivir to the WT enzyme and in the presence of mutation H274Y of H1N1 and H5N1. A long continuous MD trajectory may incur greater statistical errors as the protein denatures and evolves from one conformation to another during the time of the simulation.47 Thus, a multiple MD approach was undertaken. This method in turn reduced the force field artefacts, statistical bias, and computational time. A total simulation time of 100 ns partitioned in five distinct 20 ns MD runs was executed, with each trajectory having a different initial velocity. The multiple MD protocol applied in this study is described in Figure 3.

Figure 3

Multiple MD trajectory adopted in this report.

Abbreviation: MD, molecular dynamic.

MD simulation set-up and parameters

The MD simulation was performed using the graphics processing unit version of the PMEMD engine provided with the Amber 12 and 14 package.48,49 The FF99SB variant of the Amber force field50 was used to describe the protein. The LEAP module of Amber 12/14 allowed for the addition of hydrogen atoms to the protein and three Na+ counter ions for neutralization. The system was suspended within a TIP3P51 water box such that all protein atoms were within 8 Å of any box edge. After energy minimization, a 1,000 step restraint gradient minimization was performed. Selective boundary conditions were enforced concurrently with the particle-mesh Ewald method,52 a component of Amber 12/14, with set parameters of direct space and a van der Waals cut-off of 12 Å. This would be conducive to the treatment of long-range electrostatic interactions. The system passed through an initial energy minimization with 2,500 steps of steepest descent and a restraint harmonic potential of 500 kcal/mol Å2 being applied to the solute. After this, an additional 1,000 step unrestrained conjugate gradient energy minimization was carried out on the complete system, ie, protein, ligand, solvent molecules, and added ions. Gradual heating of the MD simulation from 0 to 300 K was executed for 50 ps, such that the system maintained a fixed number of atoms and a fixed volume, ie, a canonical (NVT) ensemble. The solutes within the system were imposed with a potential harmonic restraint of 10 kcal/mol Å2 and a collision frequency of 1.0 ps−1. Superceding the heating input, a final equilibration estimating 500 ps of the systems was required. The operating temperature was kept constant at 300 K, accompanied by the number of atoms and pressure, mimicking an isobaric-isothermal (NPT) ensemble. The systems pressure was maintained at 1 bar using the Berendsen barostat.53 The total time for each MD simulation was 20 ns, and five trajectories for each of the eight systems were generated. In each simulation, the SHAKE algorithm54 was employed to constrict the bonds of the hydrogen atoms. The time step of each simulation was 2 fs and a single-precision floating-point precision model55 was used. The simulations coincided with randomized seeding of the NPT ensemble. It was accompanied by constant pressure of 1 bar maintained by the Berendsen barostat, a pressure coupling constant of 2 ps, with a temperature of 300 K and Langevin thermostat with a collision frequency of 1.0 ps−2. Co-ordinates were saved every 1 ps and the trajectories were analyzed every 1 ps using the PTRAJ module implemented in Amber 12.0 and CPPTRAJ module in Amber 14.0.

Post-dynamic analysis

Thermodynamic calculations

An implicit solvent model was employed to describe the thermodynamic free binding energies of oseltamivir bound to WT and mutant H1N1 and H5N1 neuraminidase in this study. The molecular mechanics/generalized Born surface area method was used to evaluate the ligand-protein complex binding affinities.56–59 To calculate the free binding energy contributions, 1,000 snapshots were extracted from each of the 20 ns trajectories. The following set of equations describes the calculation of the binding free energy: where E represents the gas-phase energy, E is the internal energy, E is the Coulomb energy, and E is the van der Waals energy. The term E is directly measured from the FF99SB force field terms. The solvation energy (G) is the summation of contributions from polar and nonpolar states. The polar solvation energy contribution, G, is derived from solving the GB equation. The term G corresponds to the non-polar solvation energy contribution, which is estimated from the SASA determined using a water probe radius of 1.4 Å. The temperature and total solute entropy are represented by T and S, respectively.60

Principal component analysis

PCA reveals the structure of atomic fluctuations, and describes the motion of the system in terms of eigenvectors (planar of motion) and eigenvalues (magnitude of motion).61 The individual MD trajectories were stripped of solvent and ions using the PTRAJ and CPPTRAJ modules in Amber 12.0/14.0. The resulting trajectories were aligned against a fully minimized structure. PCA was performed on a Cα backbone with 1,000 snapshots taken every 20 ps. The first two eigenvectors (PC1 and PC2) corresponding to the first two modes of PCA covariance matrices were generated using in-house scripts. An average of the PC1 and PC2 for the 5×20 ns trajectories of the H1N1 and H5N1 WT and mutant systems was generated. The corresponding PCA scatters were plotted using Origin software (http://www.originlab.com/) and structural postscript diagrams were created using visual MDs.62 Porcupine plots of the first and second modes developed by the normal mode wizard using the ProDy interface of visual MDs were sketched for each of the systems.63

Results and discussion

MD simulations and system stability

RMSD and potential energy plots in the Supplementary materials (Figures S1–12 of H5N1 and H1N1 neuraminidase, respectively) graphically monitor the convergence of the studied systems.64 All the systems of H5N1 and H1N1 influenza viruses converge at approximately 5,000 ps by both RMSD and potential energy calculations.

RMSF and radius of gyration analysis was used to relate conformational changes and plasticity of the Cα backbone to mutation of each system.

Root mean square fluctuation

The RMSF captured the fluctuation of each residue, providing insight into the flexible regions of the protein that correspond to crystallographic β (temperature) factors.65 Figure 4 illustrates the fluctuations of systems WTH5N1 and H274YH5N1. The crystal structure of the WT neuraminidase consists of three α-helices and 24 β-strands. Evidence from the graph demonstrates an unfolding of the protein helix at residue 28 in H274YH5N1 in comparison with the WT, shortening the cardinal α-helix and β-strand by one and three residues, respectively. Closer inspection revealed a distinct change in fluctuation corresponding to region 108–144, possibly due to the presence of a hydrogen bond between Lys125 and Thr131. A similar effect was observed through β-strands 156–162 and 170–177 which share a hydrogen bond between amino acids Thr161 and Tyr171. The residue region 260–320 containing the compensatory mutation did not demonstrate large variations in fluctuations.

Figure 4

RMSF comparison of WTH5N1 and H247YH5N1: T1, T2, T3, T4, T5, and Tavg presenting the five individual 20 ns molecular dynamic trajectories and overall average, respectively.

Abbreviations: avg, average; RMSF, root mean square fluctuation; T, trajectory; WT, wild-type.

As depicted in Figure S13, by interacting with residues in the active site, the role of residue 222 appeared to shift from a structural one to a functional one. The Lys222 mutation contains a second amino group that acts as an electrophile at physiological pH; this prompts an interaction with the active site residues. Mutation I222KH5N1 caused a change in orientation of Gly332, promoting a disulfide link between parallel cysteine residues 339 and 364 and causing a shift in the β-strand orientation.66

Analysis of Figure S14 revealed a migration from a system of high flexibility (H274YH5N1) to one of low flexibility (H274Y-I222KH5N1). The change from a hydrophobic residue, His, to the aromatic Tyr274, which contains an ionizable phenolic group, introduced instability in a previously predominant hydrophobic region. Therefore, it is plausible to assume that the presence of mutation I222K in system H274Y-I222KH5N1 restricts interaction with the external solvent, thus minimizing flexibility throughout the protein.

In the H274YH1N1 system (Figure 5), an increase in fluctuation at positions 58–65, 246–260, 284, and 301–305 occurred in comparison with WTH1N1. The H274Y mutation introduced a new hydrogen bond, and the aromatic group imposed a steric bulk disrupting previous hydrophobic interactions within that region. A considerable fluctuation was observed in the latter regions of the protein. There is only a single mutation in the I222KH1N1 system (Figure S15), which induced an increase in protein flexibility because Lys enhances solvent interaction. An overlay of the RMSF plots for systems H274Y-I222KH1N1 and H274YH1N1 (Figure S16) revealed no distinct difference between the two systems, indicating that no potential compensatory effect occurred to accommodate the double mutant species of protein H274Y-I222KH1N1. Thus, the conformation of the protein remains unchanged.

Figure 5

RMSF comparison of WTH1N1 and H247YH1N1: T1, T2, T3, T4, T5, and Tavg presenting the five individual 20 ns molecular dynamic trajectories and overall average, respectively.

Abbreviations: avg, average; RMSF, root mean square fluctuation; T, trajectory; WT, wild-type.

Radius of gyration

Radius of gyration demonstrates the compactness of protein structures, providing insight into complex changes in the molecular shape.67,68 It is evident that systems WTH5N1 and H274YH5N1 (Figure 6) share a similar arrangement of amino acids in secondary and tertiary structures with an almost identical initial compactness. Although both systems share a close resemblance, the H274YH5N1 system showed a distinguishable increase in its radius of gyration over time, transmutating to a less compact structure.

Figure 6

Radius of gyration average comparison across the 20 ns molecular dynamic simulation of WTH5N1 and H247YH5N1: T1, T2, T3, T4, T5, and Tavg presenting the five individual 20 ns molecular dynamic trajectories and overall average, respectively.

Abbreviations: avg, average; T, trajectory; WT, wild-type.

System comparison of WTH5N1 and I222KH5N1 (Figure S17) unveiled the I222K mutation as the more compact protein complex. Condensation of the structure created a semipermeable complex that prevented unwarranted solvent exposure of interior amino acid residues stabilizing the mutation I222KH5N1. Comparative analysis of systems H274YH5N1 and H274Y-I222KH5N1 (Figure S18) proposed an increase in flexibility of system H274Y-I222KH5N1. The compactness of system H274Y-I222KH5N1 attempts to preserve the enzymes bioactivity whilst inferring resistance against oseltamivir. Both systems H274YH5N1 and H274Y-I222KH5N1 share a similar radius of gyration profile, which suggests that both have a similar complex compactness.

System H274YH1N1 showed a larger radius of gyration than the parent protein WTH1N1, indicating that H274YH1N1 is less tightly packed (Figure 7). Solvent interaction within the amino acid sequence perpetuated partial protein unfolding, facilitating a less compacted structure. The single mutation I222KH1N1 (Figure S19) shared a similar trend by way of possessing a large radius of gyration value as compared with the WTH1N1 system. However, isoleucine expels solvent, and when Lys was introduced, an ionizable species was generated, which, much like the Tyr residue at position 274, was capable of interacting with the surrounding solvent, relieving steric strain. The H274Y-I222KH1N1 system has a folding profile similar to that of H274YH1N1, which is indicative of an undisturbed compaction of structure. The H274Y-I222KH1N1 system (Figure S20) appeared unaffected by the mutation at position 222. This could be due to the presence of the solvent-interacting amino acid species Tyr having previously rendered the inner foldings of the protein susceptible to solvation. The introduction of an additional ionizable moiety better complemented the structure in this conformation.

Figure 7

Radius of gyration average comparison across the 20 ns molecular dynamic simulation of WTH1N1 and H247YH1N1: T1, T2, T3, T4, T5, and Tavg presenting the five individual 20 ns molecular dynamic trajectories and overall average, respectively.

Abbreviations: avg, average; T, trajectory; WT, wild-type.

Calculation of MM/GBSA free binding energy

Evaluation of the free binding energies provides a thorough insight into the total energy of each of the systems. Table 2 summarizes the average molecular mechanics/generalized Born surface area (MM/GBSA) thermodynamic energy contributions elicited from the structural data of the five 20 ns MD simulations for the H5N1 and H1N1 systems.

Table 2

Energy contribution derived from the molecular mechanics/generalized Born surface area technique corresponding to structural entities of the H5N1 system

Complexes	Trajectory	ΔGbind	ΔEele	ΔEvdW	ΔGgas	ΔGsol
H274	T1	−27.80±6.57	−54.97±15.38	−30.25±3.27	−85.21±17.19	57.42±12.03
	T2	−28.56±3.63	−57.65±7.07	−33.19±2.43	−90.85±7.05	62.29±6.22
	T3	−38.19±4.36	−71.84±6.42	−36.14±2.87	−107.97±6.83	69.79±7.34
	T4	−29.29±7.13	−60.08±8.79	−31.52±3.71	−91.63±10.82	61.34±6.48
	T5	−31.89±5.17	−62.10±10.66	−33.42±2.90	−95.52±11.96	63.64±7.78
	Tavg	−31.14±4.34	−61.33±9.16	−32.91±4.59	−94.24±8.70	62.89±7.45
I222K	T1	−24.22±6.65	−59.82±11.74	−28.20±4.20	−88.03±12.78	63.80±7.97
	T2	−23.98±5.44	−41.06±12.39	−31.92±3.75	−72.98±14.35	48.99±9.56
	T3	−23.08±4.42	−56.77±11.18	−31.16±3.60	−87.93±12.45	64.85±9.48
	T4	−26.51±7.93	−69.03±8.79	−31.44±3.22	−100.47±9.68	73.96±7.84
	T5	−22.05±8.68	−49.31±17.41	−29.29±4.90	−78.60±20.95	56.55±23.77
	Tavg	−23.97±4.89	−55.20±11.68	−30.40±6.44	−85.60±14.61	61.63±8.71
H274Y	T1	−30.80±4.62	−60.00±8.22	−31.74±2.66	−91.75±8.57	60.95±5.98
	T2	−28.53±3.97	−56.79±8.77	−29.99±3.23	−86.78±9.48	58.24±7.31
	T3	−30.80±4.62	−60.00±8.22	−31.74±2.66	−91.75±8.57	60.95±5.98
	T4	−28.53±3.97	−56.79±8.77	−29.99±3.23	−86.78±9.48	58.24±7.31
	T5	−28.53±3.97	−56.79±8.77	−29.99±3.23	−86.78±9.48	58.24±7.31
	Tavg	−29.44±3.44	−58.07±6.90	−30.69±4.11	−88.76±8.68	59.33±6.11
H274Y-I222K	T1	−16.86±4.85	−39.61±9.57	−28.15±3.75	−67.76±9.47	50.90±8.06
	T2	−19.70±5.25	−45.19±11.75	−24.53±3.98	−69.72±13.85	50.02±9.90
	T3	−15.41±6.14	−21.27±8.69	−25.64±7.34	−46.91±13.87	31.50±9.50
	T4	−20.22±3.68	−50.65±8.41	−23.39±2.77	−79.04±8.68	58.83±7.00
	T5	−19.63±3.26	−27.45±9.69	−30.82±2.68	−58.27±9.66	38.64±8.05
	Tavg	−18.36±3.98	−36.84±9.32	−26.51±5.51	−64.34±10.78	45.98±7.54

Abbreviations: ΔGbind, binding free energy; ΔEele, electrostatic energy; ΔEvdW, van der Waals energy; ΔGgas, gas-pahse energy; ΔGsol, solvation energy; avg, average; T, trajectory.

The terms ΔE and ΔE represent the electrostatic and van der Waals intermolecular interacting components, respectively, between the protein and inhibitor.69–70 Systems WTH5N1 and I222KH5N1 showed a difference of −6.1302 kcal/mol (ΔE) and −2.5055 kcal/mol (ΔE). H274YH5N1 and H274Y-I222KH5N1 showed a difference of −21.2395 kcal/mol (ΔE) and −4.1821 kcal/mol (ΔE). The energy differences advocate a more stable protein–ligand interaction of WT and H274YH5N1. This was measured in terms of effective binding distance and positioning of amino acid residues for optimum interaction between charged entities. The ΔG, difference defined between systems WTH5N1 and I222KH5N1 was estimated to 1.2626 kcal/mol and for systems H274YH5N1 and H274Y-I222KH5N1 was estimated to be 13.3466 kcal/mol. The change to polar amino acids (eg from Ile to Lys and from His to Tyr) in systems I222KH5N1 and H274Y-I222KH5N1 contributed to an improved ΔG energy difference. Closer inspection of the free binding energy, ΔG, quantified the interaction between the protein and ligand. The difference in free binding energy between systems WTH5N1 and I222KH5N1 was −7.1732 kcal/mol, and that between systems H274YH5N1 and H274Y-I222KH5N1 was −11.0751 kcal/mol. The trend in binding free energy difference was observed to be as follows: WTH5N1<I222KH5N1~H274YH5N1<H274Y-I222KH5N1, indicating that the protein-ligand interaction became progressively thermodynamically unfavorable in the presence of mutations.

Table 3 corresponds to the calculated energies for the H1N1 neuraminidase systems. The amine group of oseltamivir was kept positively charged during the MD simulation. According to the literature, binding of drug in a positively charged state to the active site residues is falsely represented by an ameliorated free binding energy when comparing the mutant sequence with the WT sequence.71 Based on a study by Le et al it was proposed that a predominantly negatively charged column of residues at the binding pocket electrostatically funnels oseltamivir to the active site of N1 neuraminidases.1 This greatly impacts the binding pose, whereby a positively charged drug is drawn toward the epicenter of the active site. The energy contribution ΔG of H274YH1N1 and I222KH1N1 is testament to the proposed binding instruction as both systems highlighted an improved binding, with an energy difference from the WT of −4.1784 kcal/mol and −2.4447 kcal/mol, respectively. Surprisingly, the double mutation species H274Y-I222KH1N1 demonstrated a similar energy profile trend to WTH1N1. The electrostatic energy of the I222KH1N1 and H274YH1N1 species differed significantly from the WTH1N1 system by −8.3000 kcal/mol and –15.9730 kcal/mol, respectively. However, a remarkable drop in electrostatic energy was observed in the H274Y-I222KH1N1 system compared with H274YH1N1. This phenomenon could relate to the conformation of the double mutant system, as the residues interacting with the solvent direct themselves inwardly, interacting with neighboring amino acid residues. The ΔE and ΔG differences between H274YH1N1 and double mutation H274Y-I222KH1N1 confirm the funneling mechanism, as both shared an improvement in the ΔE value over the WT. However, a significant decline in ΔG of 37.4933 kcal/mol indicated a lack of solvent interaction. The van der Waals contributions for I222KH1N1 suggested a slight decline in hydrophobic interaction, with an energy difference from WT of −1.0594 kcal/mol. A ΔG difference between the systems of 4.7960 kcal/mol implied enhanced solvent interaction of the I222KH1N1 complex. The H274YH1N1 offered a similar van der Waals contribution to WTH1N1, as the aromatic group replaced a linear hydrocarbon chain. Despite this structural feature, a solvation energy difference of 13.1386 kcal/mol indicated that the hydroxyl group of Tyr imposed an increased solvent exposure.

Table 3

Energy contribution derived from molecular mechanics/generalized Born surface area technique corresponding to structural entities of the H1N1 system

Complexes	Trajectory	ΔGbind	ΔEele	ΔEvdW	ΔGgas	ΔGsol
H274	T1	−20.30±4.34	−144.88±10.13	−31.31±2.52	−176.19±10.21	155.90±9.36
	T2	−16.48±9.03	−137.06±22.53	−28.17±5.26	−165.22±24.78	148.74±17.50
	T3	−25.81±4.11	−151.45±8.54	−32.68±3.16	−184.14±8.37	158.32±6.69
	T4	−30.52±4.24	−167.85±10.35	−30.61±3.50	−198.46±10.11	167.94±8.49
	T5	−23.27±4.21	−148.70±12.83	−30.35±2.82	−179.04±12.23	155.78±9.73
	Tavg	−23.28±5.19	−149.99±12.88	−30.63±3.45	−180.61±13.14	157.34±10.35
I222K	T1	−28.43±4.90	−171.09±13.93	−31.26±3.42	−202.34±13.73	173.91±11.24
	T2	−20.18±6.53	−144.76±15.31	−25.08±4.85	−169.84±16.07	149.66±12.94
	T3	−27.13±6.00	−156.67±14.22	−31.64±3.28	−188.31±15.30	161.18±11.58
	T4	−26.87±6.57	−165.80±19.73	−31.10±3.41	−196.90±20.91	170.04±16.10
	T5	−26.01±5.15	−153.12±12.22	−28.76±3.28	−181.88±13.21	155.88±10.14
	Tavg	−25.72±5.83	−158.29±15.08	−29.57±3.65	−187.85±15.84	162.13±12.40
H274Y	T1	−26.22±3.25	−131.63±16.23	−32.37±2.62	−164.00±15.71	137.78±14.03
	T2	−36.56±4.76	−175.26±15.34	−34.47±3.93	−209.72±14.01	173.16±13.30
	T3	−27.71±3.24	−179.67±11.45	−31.25±2.58	−210.91±11.38	183.21±10.88
	T4	−24.26±5.61	−173.30±16.41	−31.80±2.60	−205.10±16.83	180.84±13.83
	T5	−22.52±5.12	−169.94±22.45	−29.96±2.93	−199.91±21.39	177.39±19.32
	Tavg	−27.46±4.40	−165.96±16.37	−31.97±2.93	−197.93±15.87	170.47±14.27
H274Y-I222K	T1	−20.35±5.40	−139.80±15.52	−27.49±4.01	−167.29±15.07	146.94±12.68
	T2	−24.92±4.13	−126.91±17.15	−32.09±2.96	−158.99±15.94	134.08±13.95
	T3	−25.68±3.55	−112.27±12.68	−33.69±2.87	−145.96±11.80	120.28±10.14
	T4	−26.03±3.60	−115.52±15.64	−32.90±2.48	−148.42±15.24	122.39±13.43
	T5	−20.62±5.29	−132.49±22.51	−29.35±3.11	−161.85±22.20	141.22±19.16
	Tavg	−23.52±4.39	−125.40±16.70	−31.10±3.09	−156.50±16.05	132.99±13.88

Abbreviations: ΔGbind, binding free energy; ΔEele, electrostatic energy; ΔEvdW, van der Waals energy; ΔGgas, gas-pahse energy; ΔGsol, solvation energy; avg, average; T, trajectory.

Hydrogen bond formation between amino acid residues

The dimensionality of a protein is a vector metric governed by the direction of interaction and number of hydrogen bonds. Comparative analysis of WTH5N1 and H274YH5N1 verified a decrease in the overall number of hydrogen bonds in the mutant species (Figure 8). This phenomenon was attributed to a reduction of hydrophobic interactions within the protein, evident from the van der Waals energy contribution. The residues of neuraminidase interact to a greater extent with the surrounding solvent, hence a depletion in the number of internal protein hydrogen bonds. System I222KH5N1 showed a similar trend to the H274YH5N1 system (Figure S21). The double mutation H274Y-I222KH5N1 showed an increase in the number of hydrogen bonds (Figure S22), offering a more compact structure that is supported by the radius of gyration.

Figure 8

Averaged number of hydrogen bonds in WTH1N1 and H247YH1N1 across the 20 ns molecular dynamic simulation: T1, T2, T3, T4, T5, and Tavg presenting the five individual 20 ns molecular dynamic trajectories and overall average, respectively.

Abbreviations: avg, average; T, trajectory; WT, wild-type.

Mutation H274YH1N1 (Figure S23) showed a decline in the overall average of calculated hydrogen bonds as compared with the WTH1N1 system. This finding supports the idea of an unfolding of the protein due to depletion of interactions facilitating α-helices or β-strands capable of compacting the enzyme. Similarly, system I222KH1N1 (Figure S24), compared to the parent protein, WTH1N1, has a recorded decrease in the number of internal hydrogen bonds. A comparison of systems H274YH1N1 and H274Y-I222KH1N1 (Figure S25) revealed a near equivalent hydrogen bonding pattern. This could be due to torsional or steric stress of the macromolecule, such that thermodynamically no further unraveling of the protein can be accommodated before an irreversible deformation of the protein takes effect.

PCA is used to assess the impact of mutations on the conformational dynamics of neuraminidase to solicit an impartial resistance mechanism.72 Collation of systems WTH5N1 and H274YH5N1 (Figure 9) associated H274YH5N1 with restricted motion in relation to its α-carbon backbone. The displacements were reduced, minimizing the spatial occupancy of H274YH5N1 with a covariance of 0.3191. The magnitude of covariance for WTH5N1 was estimated to be −3.6367, advocating a more flexible framework with an inversely proportionate motion (Figure S26). System H274YH5N1 has a disproportionate motion of correlation (R2) value of 0.1814 in that the vectors commission dynamism in a single cooperative direction (Figure S27). The R2 value for WTH1N1 of −0.6526 indicated a more proportionate but antagonistic motion.

Figure 9

Principal component analysis scatter plots of 1,000 frames of the distribution along two planes, PC1 and PC2, for: WTH5N1 and H247YH5N1 illustrating differences in eigenvectors of T1, T2, T3, T4, T5, and Tavg presenting the five individual 20 ns molecular dynamic trajectories and overall average, respectively.

Abbreviations: avg, average; PC1, principal component 1; PC2, principal component 2; T, trajectory; WT, wild-type.

Examination of system I222KH5N1 (Figure S28) suggested that the neuraminidase enzyme tolerated an abridged flexibility of its α-carbon backbone in comparison with the WT enzyme. This corresponded with the determinations surmised from the RMSF values. The single mutation has a covariance of −3.1127, which is marginally smaller in magnitude than that for the WT, with an R2 of −0.3032. The direction of motion of I222KH5N1 is antiparallel, the associations predominate by an antonymous relationship of large and small variables (Figure S29). The double mutation, H274Y-I222KH5N1, has improved flexibility in the α-carbon backbone in comparison with H274YH5N1. This is demonstrated in the PCA scatter and covariance value estimated at 0.8148 (Figure S30). The improvement in flexibility coincided with the RMSD values, suggesting a greater number of conformational possibilities of the H274Y-I222KH5N1 system. An improved R2 of 0.2337, in contrast with H274YH5N1, described a unified direction of motion of disproportionate magnitude (Figure S31).

The WTH1N1 enzyme (Figure 10) appeared to have reduced internal motion of its α-carbon backbone in comparison with H274YH1N1, a restriction which corresponded to the compactness observed from calculation of the radius of gyration. The covariance and R2 values were estimated to be −1.2249 and −0.4616, respectively, suggesting a disproportionate antithetic motion of the enzyme (Figure S32). A similar effect was observed with the double mutation system H274Y-I222KH1N1, where the overall flexibility of the protein was rigid in comparison with H274YH1N1 (Figure S33). The covariance and R2 values of −0.4099 and −0.1539 support these findings (Figure S34). System H274YH1N1 demonstrated a more flexible α-carbon backbone which moved in a disorganized motion with a covariance of −0.7200 and an R2 value of −0.05218 (Figure S35). System I222KH1N1 exhibited a continuous linear motion along its α-carbon backbone, with a covariance of 0.6457 (Figure S36). This phenomenon was supported by an R2 of 0.3271 (Figure S37).

Figure 10

Principal component analysis scatter plots of 1,000 frames of the distribution along two planes, PC1 and PC2, for: WTH1N1 and H247YH1N1 illustrating differences in eigenvectors of T1, T2, T3, T4, T5, and Tavg presenting the five individual 20 ns molecular dynamic trajectories and overall average, respectively.

Abbreviations: avg, average; PC1, principal component 1; PC2, principal component 2; T, trajectory; WT, wild-type.

Solvent accessible surface area

The SASA imparts information relative to the compactness of the structure as well as the extent of hydrophobicity in the interior of the folded protein.73 System H274YH5N1 appeared to expand over time (Figure S38), whereas WTH5N1 remained unchanged. Both systems reflected a distinct and relatively stable complexed conformation, which is affirmed by the radius of gyration estimations. System I222KH5N1 (Figure S39) and H274Y-I222KH5N1 (Figure S40) demonstrated a nominally lower surface area than the WT and H274YH5N1, respectively. This implied that contact between the van der Waals partitions and the solvent was diminished. A dehydron species may be present in systems I222KH5N1 and H274Y-I222KH5N1 due to the contributions of ΔE and ΔG, resulting in a reduced interaction between the solvent and van der Waals region.74 The preservation of the hydrophobic regions results in a loss in volume of the active site.

The WTH1N1 system has a reduced exposed surface area when compared with H274YH1N1 (Figure S41). This observation corresponded to the PCA of H274YH1N1, facilitating a structure with greater potential for interaction with solvent. System I222KH1N1 exhibited a slightly different trend such that initially both the WT and single mutant shared similar SASA. Over time, at ~10,000 ps onward, a change in conformation occurred, with I222KH1N1 appearing to become more accessible to the solvent (Figure S42). A similar effect was found between H274YH1N1 and H274Y-I222KH1N1 (Figure S43). During the initial 10,000 ps of the simulation, there appeared to be a difference in the SASA of the two systems, with H274YH1N1 having greater solvent accessibility than the double mutant. In the later 10,000 ps, both systems begin to overlap. This anomaly suggests a conformational change in the presence of the double mutation H274Y and I222K, whereby these residues pointed outward, encouraging solvent interaction and disrupting the three-dimensional compactness.

Conclusion

Significant findings emerged in our endeavor to gain insight into the binding mode and origin of resistance of oseltamivir in H1N1 and H5N1. It was noted that the mechanism of resistance is entirely dependent on the type of mutation, relative positioning of such occurrences, and the repercussions thereof on the active site of the enzyme under analysis. Our results indicate that resistance of neuraminidase was inherited by the expression of H274Y and/or I222K mutations. The neuraminidase of H1N1 has known resistance against oseltamivir in the presence of mutations H274Y and/or I222K. Despite originating from the same host by reassortment, H5N1 has proven to be more susceptible to oseltamivir than H1N1 in the presence of mutations. Comparatively, neuraminidase of H5N1 appears to experience greater structural fluctuations in acquiring resistance against oseltamivir as opposed to the H1N1 system. H5N1 gains resistance by a loss in volume of the active site due to compaction of the neuraminidase, while H1N1 develops resistance by solvent exposure of active site residues. Subsidiary analyses of RMSD, RMSF and potential energy confirm the progressive development of resistance from WT through to mutants (for both H1N1 and H5N1).

The free binding energy established from the MM/GBSA algorithm offers the first sign of evidence of resistance. Based on the electrostatic funnel mechanism of neuraminidase (for H1N1), an increase in binding energy was observed. This was an artefact of the active site residues and drug interacting with solvent, rather than with one another. The presence of the H274Y mutation has a compensatory effect to ensure the survival of the viral species by inhibiting enzyme interaction with the drug, without compromising the integrity of the enzyme. Mutation I222K, through distance, transposes a putative constriction in the volume of the active site of H5N1, whereas in H1N1 the protein unfolds, promoting solvent interaction. As a result, the active site inherently cannot support the bulky 1-ethylpropoxy hydrophobic moiety of oseltamivir. It is reasonable to suppose that the I222K mutation has a potent effect in discriminating between binding of the substrate and binding of the drug oseltamivir. The H274Y-I222K double mutant exhibits even greater drug resistance. Overall resistance in systems H5N1 and H1N1 occurs as a result of loss of hydrophobicity within the binding pocket of the active site. The lack of hydrophobicity leads to structural collapse of the available active site, which in turn diminishes the integrity of the binding landscape. This poses a threat, given that H5N1/avian influenza has been found to have an increasing incidence of human transmissibility. The existence of a resistant strain is having a severe negative impact on human health. Therefore, it is imperative to design and identify unique, innovative, and selective chemotherapeutic agents to adapt to the newly defined binding pocket.