Single H5N1 influenza A neuraminidase mutation develops resistance to oseltamivir due to distorted conformational and drug binding landscape: multiple molecular dynamics analyses

Abstract

Clinical studies showed that a single mutation, I117V, develops severe resistance to oseltamivir, the first orally active influenza A neuraminidase inhibitor, in highly pathogenic H5N1 influenza A viruses. However the molecular mechanism by which this single mutation confers oseltamivir resistance is not well understood. Herein, and via “multiple” molecular dynamics approaches, we conducted a comparative study on the wild type and I117V mutant to understand the origin of resistance at the molecular level using a wide range of post-dynamics analytical tools. These include binding free energy calculations (MM/GBSA), principal component analysis (PCA), radius of gyration (RG), per-residue fluctuation (RMSF), solvent accessible surface area (SASA) and per-residue secondary structure occupancy. Results showed that single mutations at position 117 of H5N1 caused (1) distortion of the orientation of oseltamivir in the active site, (2) decrease in the binding affinity by ∼6 kcal mol⁻¹ when compared to the wild type, (3) variation in the overall enzyme motion as evident from the PCA for both systems and (4) distortion of the hydrogen bonding network and atomic interactions with the inhibitor.

---

1. Introduction

The highly pathogenic influenza A (H5N1) predominantly originates from avian species but intermittently crosses the species barrier to infect humans to cause drastic respiratory complications and death. Cases of human H5N1 virus infection were first reported in 1997 and the virus continues to cause sporadic human infections worldwide with a total 649 reported cases and 385 reported deaths.¹ Regardless of the incapability of the avian H5N1 virus to transmit from human to human, potential emergence of a pandemic virus is still an inevitable threat. Continuous evolution and genetic diversity of the circulating H5N1 virus is a major concern as it poses a potential threat to human health and could result in the emergence of pandemic human H5N1 strains.²

Influenza viruses launch cellular invasion through two significant glycoproteins, haemagglutinins and neuraminidase. Haemagglutinins intercedes cell-surface sialic acid receptor adhesion to initiate virus infection.³ Post viral replication, neuraminidase subsequently cleaves sialic acid from the nascent virion thus permitting the release and spreading of influenza virus to uninfected cells.³ This significant role of neuraminidase led to its identification as the drug target.⁴ Neuraminidase inhibitors arrest the hydrolysis of sialic acid linkages, particularly the adjacent and N-acetylated carbohydrate molecule of cellular glycoprotein.⁵

Human infections by H5N1 influenza are predominantly caused by clades 0, 1, 2.1, 2.2, 2.3 and 7 since 1997 but 2010–2012 infections were from clades 2.3.2 and 2.3.4.^2,6 Efforts to develop an effective cross-clade H5N1 immunogenic vaccine in response to this threat is pending, meanwhile antiviral prophylaxis and treatment plays a pivotal role in the absence of a potent vaccine. Oseltamivir (Tamiflu) and zanamivir (Relenza) are the current neuraminidase inhibitors for treatment of influenza A and B infections.⁷ Oseltamivir preference as the main antiviral drug for influenza treatment and prevention is due to its oral bioavailability and has been proven active in vitro against human and avian H5N1 viruses, however isolates housing neuraminidase point mutations may reduce oseltamivir sensitivity.⁸ Oseltamivir is a transition state sialic acid analogue with a bulky hydrophobic group, which is lacking in sialic acid, designed to target influenza N2 neuraminidase. The bulky side chain forces the active site amino acid rearrangements to accommodate oseltamivir binding.⁹ Alarming resistance to this drug was initially attributed to the global spreading of oseltamivir-resistant H1N1 influenza strains harboring mutations conferring resistance.^10,11 Some of these NAI resistant mutations have been identified in HPAI H5N1 virus isolates as well, including H274Y and N294S N1 neuraminidase subtype mutations as well as E119V, R292K and N294S in N2 subtype which confer oseltamivir resistance.¹² Oseltamivir resistant clones harbouring an I117V mutation were also detected from the human H5N1 virus in Turkey in 2009,¹ Bangladesh and India in 2010–2011.¹³ This mutation was determined to slightly reduce oseltamivir susceptibility in vitro and dramatically in vivo, which was attributed to the hydrogen bond loss between R118 and oseltamivir carboxylate thus resulting in reduced binding affinity.¹³

This is of global concern as the entrusted neuraminidase inhibitors could become completely ineffective as primary influenza treatment in future, should virulent resistant influenza pandemic strains emerge in the absence of a vaccine. Previous studies have attributed oseltamivir resistance to hydrophobic interactions of lipophilic side-chains with influenza A and B neuraminidase with homologous active site. H5N1 resistance to oseltamivir has been defined to be caused by mutations at residue positions 119, 152, 274 and 292.^14,15 It was also proposed that the active site undergoes reorientation to accommodate the bulky side chain of oseltamivir.¹⁶ However, such proposal was developed from the X-ray crystal structure of N9 neuraminidase–oseltamivir complex.¹⁶ Therefore, thorough understanding of the mechanism of resistance of such variants is crucial to the development of effective and potent antiviral drugs. Although the properties of influenza H5N1, neuraminidase subtype single and dual mutations have been extensively studied as well as the examination of the viral strength against neuraminidase inhibitors,^13,14,17,18 however, features of the I117V mutant (Fig. 1) and its contribution to oseltamivir resistance has not been deeply examined.

Fig. 1 The structure of H5N1 influenza A virus neuraminidase showing the position of the mutation (red), I117V, studied in this work.

Attempts to obtain molecular understanding on how this mutation confers resistance to oseltamivir have been previously conducted.^1,10,19 However, findings from such studies result in postulations of drug resistance mechanisms that require prolonged timescale MD of neuraminidase structure but are imperatively restricted by the transient time scales that can be fulfilled by these methods.²⁰ Advanced MD calculations and analysis have not been conducted on H5N1 neuraminidase–oseltamivir I117V mutant. In this study, we aim to provide deeper insight into the effect of this neuraminidase mutation and adding more dimensions to the experimental work that has been previously conducted. To achieve this, we performed advanced molecular dynamics simulations of oseltamivir–neuraminidase I117V mutant as well as wild-type neuraminidase, to further examine molecular dynamics affecting oseltamivir binding. Findings from this study could be valuable in the future development of more vigorous neuraminidase inhibitors.

Molecular dynamics simulations presents conformational changes of a single molecule which is examined over time, whereas in the actual biological system, molecule conformational changes are demonstrated collectively.²¹ This is attributed to the emergence of different pathways displaying similar free-energy barriers in the unfolding of a protein.^22,23 Therefore, applying multiple-trajectory molecular dynamics simulations could sample conformational space more efficiently by minimizing the impact of entrapment in local minima, thus providing more relevant and reliable insight to protein dynamics.^21,24

Diverse post-dynamics techniques have been applied to obtain molecular understanding of molecular dynamics simulations. Principal component analysis (PCA), alternatively known as essential dynamics analysis, is one of the widely applied post-dynamics techniques to comprehend the behavior in biological systems.²⁵ PCA eliminates comprehensive translational and rotational motions in the MD trajectory and presents concerted motions necessary to understand.^26–28 This technique describes collective atomic displacements and can identify major conformational changes between the structures and has been previously employed to pinpoint the difference in motions of the wild-type and mutant strains.^25,29 In this study we perform 100 ns multiple MD simulations, PCA and binding free energy calculations to gain insight into the effect of I117V mutation on oseltamivir resistance. The coherent computational methods applied in this study could shine light on the effect of mutations on drug resistance and also assist in drug development and discovery of potent antiviral drugs.

2. Computational methods

2.1. Systems preparation

The X-ray crystal structure of wild-type influenza N1 neuraminidase, PDB code 2HU4 ³ in complex with oseltamivir was obtained from Protein Data Bank. Neuraminidase X-ray crystal structure exists as a homo-tetramer but only one chain (chain A) was used for simulations in this study to reduce computational cost. A mutation of isoleucine to valine at position 117 (I117V) was performed in silico from the wild-type crystal structure using Chimera software package.³⁰ Ligand and receptor modification as well as visualizations were conducted in Chimera and Avogadro software³¹ respectively.

2.2. Molecular dynamics simulations

To examine the effect of a single mutation I117V on oseltamivir binding with neuraminidase, a multiple-MD trajectory approach was performed. Since a continuous MD trajectory approach could result in discrepancies such as high statistical errors arising from denatured protein configuration during a simulation, a multiple-MD trajectory approach is therefore more reliable with reduced force field-induced discrepancies, statistical bias and computational duration. In this study, a multiple-MD trajectory approach was performed from an initial configuration generated from a 2 ns MD run, followed by 4 distinctive MD runs for a period of 25 ns with different respective initial velocities as described in Fig. 2.

Fig. 2 Graphical summary of multiple-MD trajectory approach applied in this study.

2.2.1 Molecular dynamics simulation set-up and parameters. Simulations of neuraminidase–oseltamivir complexes (wild-type and mutant) were performed on the GPU version of PMEMD packaged in Amber12.^32,33 The AMBER force field³⁴ FF99SB variant was applied to describe the protein. The LEAP module implemented in Amber12 was used to add hydrogen atoms to the protein and to add counter ions for systems neutralization. The respective systems were contained in a TIP3P water box³⁵ with protein atoms located 8 Å away from the water box edge. The periodic boundary conditions were adopted and long-range electrostatics were treated with particle-mesh Ewald method³⁶ implemented in Amber12 with direct space and a 12 Å van der Waals cut-off. Initial energy minimization of the respective systems was performed with as restraint potential of 500 kcal mol⁻¹ Å⁻² applied to the solute, for 1000 steps of steepest descent followed by conjugate gradient minimization of 1000 steps. Unrestrained conjugate gradient minimization of 1000 steps was performed for the entire system. Canonical ensemble (NVT) MD simulations were performed for 50 ps with gradual heating from 0 to 300 K, harmonic restraints of 5 kcal mol⁻¹ Å⁻² for solute atoms and a Langevin thermostat with a 1 ps random collision frequency. The systems were equilibrated at 300 K in the NPT ensemble for 500 ps without imposed restraint. A Berendsen barostat was used to maintain the pressure at 1 bar. All hydrogen bonds were constrained using the SHAKE algorithm,³⁷ and a timescale of 2 fs for all MD runs using the SPFP precision model³⁸ was applied. A short production run of 2 ns was performed prior to a multiple MD production run of 25 ns × 4 from configurations with Boltzmann weighted random velocities extracted at time intervals of 500 ps, 1 ns, 1.5 ns and 2 ns from the initial 2 ns MD production run with an isothermal isobaric (NPT) ensemble using a Berendsen barostat³⁹ for each case with a target pressure of 1 bar and a 2 ps pressure coupling constant. The coordinates were saved every 1 ps time interval and the trajectories were analyzed every 1 ps and further analyses including RMSF, RMSD, radius of gyration, solvent accessible surface area, hydrogen bond occupancy, secondary structure analysis and PCA were performed using the Amber12-implemented modules, PTRAJ and CPPTRAJ.⁴⁰ Visualizations and plots were conducted using Chimera molecular modeling tool for data analysis, respectively.

2.3. Thermodynamic calculations

The binding free energies of oseltamivir bound neuraminidase wild-type and mutants were calculated using the Molecular Mechanics/Generalized-Born Surface Area method.^41–46 Binding free energy calculation is an endpoint energy calculation which provides valuable information about ligand–complex association. Calculation of binding free energies considered 1000 snapshots from each 25 ns trajectory. Average values of the 4 trajectories were computed (see Table 2). The following set of equations describes the binding free energy calculation:

ΔG_bind = G_complex − G_receptor − G_ligand (1)

ΔG_bind = E_gas + G_sol − TS (2)

E_gas = E_int + E_vdw + E_ele (3)

G_sol = G_GB + G_SA (4)

G_SA = γSASA (5)

where E_gas is the gas-phase energy; E_int is the internal energy; E_ele are the Coulomb and E_vdw are the van der Waals energies. E_gas is directly estimated from the FF99SB force field terms. The salvation free energy, G_sol, can be decomposed into polar and non-polar states of contribution. The polar salvation, G_GB, contribution is estimated by solving the GB equation whereas the non-polar salvation contribution, G_SA, is determined from the solvent accessible surface area (SASA) estimated by using a water probe radius of 1.4 Å. T and S are the temperature and the total solute entropy, respectively. To evaluate the individual amino acid contribution to the total binding free energy between the inhibitor oseltamivir and neuraminidase (wild-type and mutant), the interaction energy decomposition analysis per residue was computed using MM/GBSA binding free energy method in Amber12.

Table 1 Average hydrogen bond distances between residue atoms interacting with oseltamivir. Atom labels are provided in the H-bonding interaction network plot in the ESI material (S6)

H-bond	Wild	Mutant
	Distance (Å)	Distance (Å)
Glu37 (OE2)⋯(N2) Oselt	4.01	—
Asp69 (OD1)⋯(N2) Oselt	5.07	5.31
Arg70 (NH2)⋯(O3) Oselt	3.46	3.75
Arg211 (NH1)⋯(O1) Oselt	2.91	3.40
Arg211 (NH2)⋯(O1) Oselt	3.06	3.44
Tyr262 (OH)⋯(O1) Oselt	3.48	4.77
Arg286 (NH1)⋯(O1) Oselt	2.98	3.55
Arg286 (NH2)⋯(O2) Oselt	3.09	3.65

---

2.4. Principal component analysis

Prior to MD trajectory processing for PCA, individual 25 ns MD trajectories of all complexes were stripped of the solvent and ions using the PTRAJ⁴⁷ module in Amber12. The stripped trajectories were subsequently aligned against a fully minimized structure. PCA was performed on C-α atoms on 1000 snapshots at 5 ps time interval each. The first two principal components (PC1 and PC2) were computed and the covariance matrices were generated using in-house scripts. PC1 and PC2 generated from the 4 individual 25 ns trajectories were averaged for both wild-type and mutants. PC1 and PC2 correspond to the covariance matrix first two Eigenvectors. PCA scatter plots were then constructed in Origin software (http://www.originlab.com/) and the structural diagrams were constructed in VMD.⁴⁸ Porcupine plots were constructed using ProDy interface of VMD Normal Mode Wizard (NMW).⁴⁹

3. Results and discussion

3.1. Systems stability and MD simulations

To ascertain the equilibration of the systems prior to MD analysis, the root of mean square deviation (RMSD) and potential energy fluctuation were monitored. Plots of RMSD and potential energy as well as PDB coordinates of the simulated systems are provided in the ESI section.†

3.2. Post-dynamics analysis: wild-type versus I117V mutant

3.2.1 Root of mean square fluctuation (RMSF). Fig. 3 presents root of mean square fluctuations of amino acid residues from wild-type and I117V mutant simulations. An intriguing observation is that the amino acid residues at region 100–150 (a region possessing the I117V mutation), shows higher mutant fluctuation as compared to the wild-type. This is explained by an apparent observation that the valine residue at region 117 is one methylene-group smaller than the substituted isoleucine and thus induces a slight distortion around neighbouring residue 117 regions. Although the overall structures of wild and mutant show a relatively similar trend of amino acid fluctuations, however, relatively higher amino acid fluctuations are observed in the I117V mutant neuraminidase structure in comparison to the wild-type. This phenomenon could have an influence the resulting binding affinity.

Fig. 3 RMSF plots of the wild-type and I117V mutant neuraminidase–oseltamivir complex systems: T1, T2, T3, T4 and T_avg presents the 4 individual 25 ns MD trajectories and the overall average, respectively. A zoomed-out view of the region containing the point mutation (100–150) is shown in T_avg plot.

3.2.2 Radius of gyration. The radius of gyration is a parameter linked to the tertiary structural volume of a protein and has been applied to obtain insight into the stability of molecules in a biological system along the MD simulation.⁵⁰ Here we evaluated the radius of gyration for the wild and mutant neuraminidase structures. As evident in Fig. 4, the radius of gyration for the I117V mutant is significantly higher in comparison to the wild-type throughout the simulation period. The radius of gyration for the mutant structure increased drastically from the beginning of the simulation up to 20 000 ps, implying a dramatic expansion in protein structure up to this point of time (20 000 ps), thus creating ample vacant room for solvent molecules to penetrate the hydrophobic domain. This imposes increased hydrophobicity to the protein surfaces thus inducing the ease of residue fluctuation around the active site and eventually resulting in reduced oseltamivir binding affinity.

Fig. 4 Radius of gyration plots of the wild-type and I117V mutant neuraminidase–oseltamivir complex systems: T1, T2, T3, T4 and T_avg presents the 4 individual 25 ns MD trajectories and the overall average, respectively.

3.2.3 Solvent accessible surface area (SASA). It was also evident in Fig. 5 that the SASA for the mutant structure is higher compared to that of the wild-type. Mutant structure SASA shows a dramatic increase along the trajectory and slightly levels after 15 000 ps as compared to the wild-type, thus implying that the structure further unfolds throughout the simulation period. Such event results in the loss of structural integrity and activity therefore causing the drug to bind inappropriately in the active site. Therefore the SASA induced by the I117V mutation markedly influences the structure as well as activity of neuraminidase.

Fig. 5 SASA of wild-type and mutant NA–oseltamivir complex systems: T1, T2, T3, T4 and T_avg present the 4 individual 25 ns trajectories and overall average, respectively.

3.2.4 Hydrogen bond formation between amino acid residues. Inter-amino acid residue hydrogen bond formation is a major steering force in controlling protein structure conformation. The introduction of a mutation in a protein is expected to result in hydrogen bond changes around the site of mutation. Therefore we examined hydrogen bond formation during the course of a simulation for both the oseltamivir-bound wild-type neuraminidase and I117V mutant neuraminidase complexes (Fig. 6). The I117V mutant complex displayed a relatively lower hydrogen bond formation during the simulation compared to the wild-type complex. The loss of hydrogen bonds leads to structural instability and affects its conformation and eventually drug binding. It was also evident that the loss in hydrogen bonding also contributed to the conformational distortion of oseltamivir and drug binding landscape (see S1†).

Fig. 6 Number of hydrogen bond formation during a simulation over time between oseltamivir-bound wild-type and I117V mutant neuraminidase complex.

We also examined the interaction of amino acid residues with oseltamivir in the active site of wild-type and I117V mutant neuraminidase. These interactions are presented in S6 in ESI material,† where greater amino acid interactions with oseltamivir are seen in the wild-type compared to the I117V mutant, which doesn't interact with Glu37. This implies a loss of a hydrogen bond between oseltamivir and Glu37 thus resulting in oseltamivir distortion and reduction of the binding affinity of oseltamivir to the active site.

Hydrogen bond network. To investigate the impact of the I117V mutation on oseltamivir binding, we monitored the evolution of hydrogen bond distances between amino acid residues interacting with oseltamivir in the active site of the wild-type and I117V mutant neuraminidase complexes throughout the simulation (Table 1). An increase in distance between hydrogen bonds of interacting residue atoms and oseltamivir was observed in a mutant complex, with ∼0.55 average difference compared to wild-type. This difference in H-bond distance is in accordance with the estimated oseltamivir binding affinity. Detailed plots of the results are presented in the ESI material section (see S4 and S5†). A summary of the average distances obtained throughout the simulation period is presented in Table 1.

3.2.5 Per-residue secondary structure elements. Per-residue occupancy of secondary structure elements were studied to obtain more insight on the structural plasticity of the respective wild-type and the I117V mutant neuraminidase structures. As evident in Fig. 7 the average occupancy of secondary structural components such as alpha-helix, 3₁₀-helix, parallel beta-sheet, antiparallel beta-sheet and turn were reduced in the I117V mutant complex structure in comparison to the wild-type complex however, a minor difference was observed between the two complexes in the case of antiparallel beta-sheet component. Pi-helix was the only secondary structure element of the I117V mutant complex that was dominant over the wild-type complex. Secondary structural elements per-residue occupancy further noted the differences between the wild and I117V mutant secondary structural elements.

Fig. 7 Comparative per-residue secondary structure occupancy between the native and I117V variant neuraminidase complex. (A–F) denotes parallel beta-sheet, antiparallel beta-sheet, 3₁₀-helix, alpha-helix, pi-helix and turn respectively.

3.2.6 MMPB-SA binding free energy calculation. Molecular mechanics and solvation energy components were averaged over respective 4 × 25 ns multiple MD simulations by applying MM/GBSA technique as listed in Table 1. The estimated binding free energy between oseltamivir and wild-type is −33.4533 kcal mol⁻¹ compared to −26.6789 kcal mol⁻¹ of the I117V mutant. This reduction in binding affinity (−6.7744 kcal mol⁻¹) due to mutation could markedly impair drug binding and reduce drug binding affinity as well as the efficiency of the drug towards the mutant. The reduction in binding energy is also attributed to the loss of the bond between Glu37 OE2 atom and oseltamivir N₂ atom (see S6†) in the active site of I117V neuraminidase mutant. These results are congruent to the experimental data – (IC₅₀ value for the I117V mutant is ∼5 times lower than the wild-type).⁵¹ Contribution of the calculated van der Waals (ΔE_vdW) to the binding free energy of the oseltamivir-bound wild-type neuraminidase complex (−33.4533 kcal mol⁻¹) are evidently higher than that of the oseltamivir-bound I117V mutant neuraminidase complex (−26.6789 kcal mol⁻¹). The estimated electrostatic contributions (ΔE_ele) to the binding free energy of the oseltamivir-bound I117V mutant neuraminidase complex (−60.2881 kcal mol⁻¹) are lower compared to those of the oseltamivir-bound wild-type neuraminidase complex (−70.9216 kcal mol⁻¹). In the oseltamivir-bound wild-type neuraminidase complex, the calculated solvation energy contributions (ΔG_sol) to the binding free energy (70.6804 kcal mol⁻¹) are higher than those of the mutant oseltamivir-bound I117V mutant neuraminidase complex (63.4584 kcal mol⁻¹). The dissected free energy components presented in Table 2 are suggestive of the ΔE_vdW and ΔE_ele to be the major energy contributors to oseltamivir binding. A significant difference in ΔG_gas (13.9963 kcal mol⁻¹) due to the I117V mutation implies that oseltamivir could lose its efficacy due to this mutation.

Table 2 MM/GBSA binding free energies of oseltamivir bound to wild-type, I117V mutant on H5N1 neuraminidase

Complexes		ΔGbind	ΔEele	ΔEvdW	ΔGgas	ΔGsol
WT	T1	−28.6905 ± 0.3726	−64.1037 ± 0.9341	−29.8198 ± 0.2879	−93.9235 ± 1.0927	65.2330 ± 0.8537
	T2	−34.0956 ± 0.2703	−74.6245 ± 0.6034	−35.1750 ± 0.2610	−109.7995 ± 0.6411	75.7039 ± 0.4983
	T3	−36.1516 ± 0.4404	−74.9884 ± 0.5659	−34.4623 ± 0.2899	−109.4507 ± 0.6484	73.2991 ± 0.3978
	T4	−34.8756 ± 0.4228	−69.9699 ± 0.7292	−33.3915 ± 0.2517	−103.3613 ± 0.8019	68.4857 ± 0.5001
	Tavg	−33.4533	−70.9216	−33.2121	−104.1337	70.6804
I117V	T1	−29.1745 ± 0.3546	−65.4665 ± 0.8735	−31.0045 ± 0.2590	−96.4709 ± 0.9445	67.2964 ± 0.7254
	T2	−23.7009 ± 0.7817	−53.8877 ± 1.1827	−29.0980 ± 0.5100	−82.9857 ± 1.5601	59.2847 ± 0.9006
	T3	−28.9884 ± 0.5961	−61.8464 ± 1.2755	−33.1625 ± 0.3123	−94.7109 ± 1.4480	65.7225 ± 0.9232
	T4	−24.8519 ± 0.6217	−59.9520 ± 1.0601	−26.4303 ± 0.5491	−86.3822 ± 1.4479	61.5303 ± 1.0223
	Tavg	−26.6789	−60.2881	−29.9237	−90.1374	63.4584

---

3.2.7 Per-residue contribution to binding free energies. We further decomposed the binding free energy into contributions from individual amino acid residues of neuraminidase. Table 2 presents protein–ligand interactions between the wild-type and I117V mutant. It is evident from the energy decomposition analysis that in the oseltamivir-bound wild-type neuraminidase system, major contributions were −11.35, −7.41 and −7.05 kcal mol⁻¹ from amino acid residues Arg286, Arg211 and Arg70 respectively. As evident in Table 2, the decomposed van der Waals and electrostatic energies in the oseltamivir-bound I117V mutant neuraminidase complex, from amino acid residues Val117, Arg286, Arg70, Arg211 and Tyr320 implies that the change in the van der Waals is partly responsible for interaction energy depreciation. Moreover, when isoleucine mutates to valine at position 117, the van der Waals contribution to the total binding energy from this residue slightly increased by one-fold whereas the electrostatic contribution decreased by 7-fold. Therefore, the mutation results in increased van der Waals and decreased electrostatic interactions (Table 3) thus affecting the total binding energy.

Table 3 Decomposed van der Waals and electrostatic energies (kcal mol⁻¹) for residues 70, 211, 286, 320 and 117 in oseltamivir-bound wild-type and I117 mutant complexes (values were averaged over the 4 individual trajectories)

Residues	van-der Waals	Electrostatic
a Mutant.
Arg70	−0.72 ± 0.45^a	−5.68 ± 4.32^a
	−0.78 ± 0.58	−7.05 ± 1.70
Val117 I117	−0.006 ± 0.00^a	0.059 ± 0.02^a
	−0.005 ± 0.00	0.066 ± 0.02
Arg211	−1.09 ± 0.76^a	−6.17 ± 2.39^a
	−1.17 ± 0.51	−7.41 ± 1.28
Arg286	−0.03 ± 0.68^a	−10.14 ± 2.37^a
	−0.05 ± 0.51	−11.35 ± 2.02
Tyr320	−1.29 ± 0.31^a	−0.69 ± 0.41^a
	−1.41 ± 0.26	−0.73 ± 0.43

---

3.2.8 Principal component analysis (PCA). The flexibility of oseltamivir-bound wild-type NA and oseltamivir-bound I117V mutant NA complexes was assessed using PCA method. Fig. 8 shows a scatter plot constructed for the wild-type and I117V mutant depicting a significant difference between the two systems plotted along the path of two principal components. As evident in Fig. 8, the I117V mutant complex occupies a larger phase space and exhibits a higher fluctuation compared to the wild-type complex. Such evidence implies that the mutant structure had a high degree of flexibility such that it was incapable to bind with oseltamivir.

Fig. 8 PCA scatter plot of 1000 snapshots along the pair of first two principal components, PC1 and PC2 for wild-type and I117V mutant showing difference in eigenvectors. Eigen values were averaged over the 4 individual 25 ns MD trajectories.

4. Conclusion

A detailed molecular understanding of the impact of a single mutation I117V on highly pathogenic H5N1 influenza neuraminidase is lacking. In this study we explored numerous computational approaches with an intention to provide a multidimensional view on the resistance of I117V mutation towards oseltamivir. These approaches, including multiple MD simulations, principal component analysis, binding free energy calculations, radius of gyration and solvent accessible surface area, assisted us to understand the impact of the I117V mutation on drug resistance. Our findings show that this mutation reduced oseltamivir binding affinity by ∼6 kcal mol⁻¹; there's a loss of overall intermolecular hydrogen bonds in the mutant complex resulting to a loose structural conformation; residue-based fluctuations are affirmative of the distorting effect of the I117V mutation on oseltamivir–neuraminidase complex; a large radius of gyration of the mutant complex decreases the interaction among neighbouring amino acid residues resulting in an unstable moment of inertia; a large solvent accessible surface area of the mutant complex permits solvent molecules to invade the hydrophobic region thus interrupting the drug binding process. This proves the instability of the mutant complex and is justified by RMSF and radius of gyration findings.

Conflict of interest

Authors declare no potential financial and other conflict of interests.

Acknowledgements

The authors acknowledge the UKZN School of Health Sciences, NRF for funding and the Center for High Performance Computing (CHPC, http://www.chpc.ac.za) for computational resources.

References

1. E. A. Govorkova, N. A. Ilyushina, J. L. McClaren, T. S. P. Naipospos, B. Douangngeun and R. G. Webster, Susceptibility of Highly Pathogenic H5N1 Influenza Viruses to the Neuraminidase Inhibitor Oseltamivir Differs In Vitro and in a Mouse Model, Antimicrob. Agents Chemother., 2009, 53, 3088–3096.
2. E. A. Govorkova, T. Baranovich, P. Seiler, J. Armstrong, A. Burnham, Y. Guan, M. Peiris, R. J. Webby and R. G. Webster, Antiviral resistance among highly pathogenic influenza A (H5N1) viruses isolated worldwide in 2002–2012 shows need for continued monitoring, Antiviral Res., 2013, 98, 297–304.
3. R. J. Russell, L. F. Haire, D. J. Stevens, P. J. Collins, Y. P. Lin, G. M. Blackburn, A. J. Hay, S. J. Gamblin and J. J. Skehel, The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design, Nature, 2006, 443, 45–49.
4. S. Q. Wang, Q. S. Du and K. C. Chou, Study of drug resistance of chicken influenza A virus (H5N1) from homology-modeled 3D structures of neuraminidases, Biochem. Biophys. Res. Commun., 2007, 354, 634–640.
5. D. J. Bucher and E. Kilbourn, A2 (N2) Neuraminidase of X-7 Influenza-virus recombinant – Determination of molecular-size and subunit composition of active unit, J. Virol., 1972, 10, 60–66.
6. J. L. McKimm-Breschkin, P. W. Selleck, T. B. Usman and M. A. Johnson, Reduced sensitivity of influenza A (H5N1) to oseltamivir, Emerging Infect. Dis., 2007, 13, 1354–1357.
7. L. V. Gubareva, L. Kaiser and F. G. Hayden, Influenza virus neuraminidase inhibitors, Lancet, 2000, 355, 827–835.
8. M. D. de Jong, T. T. Thanh, T. H. Khanh, V. M. Hien, G. J. D. Smith, N. V. Chau, B. V. Cam, P. T. Qui, D. Q. Ha, Y. Guan, J. S. M. Peiris, T. T. Hien and J. Farrar, Brief report – Oseltamivir resistance during treatment of influenza A (H5N1) infection, N. Engl. J. Med., 2005, 353, 2667–2672.
9. C. U. Kim, W. Lew, M. A. Williams, H. T. Liu, L. J. Zhang, S. Swaminathan, N. Bischofberger, M. S. Chen, D. B. Mendel, C. Y. Tai, W. G. Laver and R. C. Stevens, 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.
10. Y. Matsuzaki, K. Mizuta, Y. Aoki, A. Suto, C. Abiko, K. Sanjoh, K. Sugawara, E. Takashita, T. Itagaki, Y. Katsushima, M. Ujike, M. Obuchi, T. Odagiri and M. Tashiro, A two-year survey of the oseltamivir-resistant influenza A(H1N1) virus in Yamagata, Japan and the clinical effectiveness of oseltamivir and zanamivir, Virol. J., 2010, 7, 53.
11. N. J. Dharan, L. V. Gubareva, J. J. Meyer, M. Okomo-Adhiambo, R. C. McClinton, S. A. Marshall, K. S. George, S. Epperson, L. Brammer, A. I. Klimov, J. S. Bresee and A. M. Fry, Infections With Oseltamivir-Resistant Influenza A(H1N1) Virus in the United States, JAMA, J. Am. Med. Assoc., 2009, 301, 1034–1041.
12. J. L. McKimm-Breschkin, Influenza neuraminidase inhibitors: antiviral action and mechanisms of resistance, Influenza Other Respir. Viruses, 2013, 7, 25–36.
13. R. Takano, M. Kiso, M. Igarashi, L. Quynh Mai, M. Sekijima, K. Ito, A. Takada and Y. Kawaoka, Molecular Mechanisms Underlying Oseltamivir Resistance Mediated by an I117V Substitution in the Neuraminidase of Subtype H5N1 Avian Influenza A Viruses, J. Infect. Dis., 2013, 207, 89–97.
14. M. Z. Wang, C. Y. Tai and D. B. Mendel, 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.
15. M. A. Rameix-Welti, F. Agou, P. Buchy, S. Mardy, J. T. Aubin, M. Veron, S. van der Werf and N. Naffakh, Natural variation can significantly alter the sensitivity of influenza A (H5N1) viruses to oseltamivir, Antimicrob. Agents Chemother., 2006, 50, 3809–3815.
16. A. Moscona, Drug therapy – Neuraminidase inhibitors for influenza, N. Engl. J. Med., 2005, 353, 1363–1373.
17. N. Renzette, D. R. Caffrey, K. B. Zeldovich, P. Liu, G. R. Gallagher, D. Aiello, A. J. Porter, E. A. Kurt-Jones, D. N. Bolon, Y.-P. Poh, J. D. Jensen, C. A. Schiffer, T. F. Kowalik, R. W. Finberg and J. P. Wang, Evolution of the Influenza A Virus Genome during Development of Oseltamivir Resistance In Vitro, J. Virol., 2014, 88, 272–281.
18. M. L. Mihajlovic and P. M. Mitrasinovic, Another look at the molecular mechanism of the resistance of H5N1 influenza A virus neuraminidase (NA) to oseltamivir (OTV), Biophys. Chem., 2008, 136, 152–158.
19. A. C. Hurt, P. Selleck, N. Komadina, R. Shaw, L. Brown and I. G. Barr, Susceptibility of highly pathogenic A(H5N1) avian influenza viruses to the neuraminidase inhibitors and adamantanes, Antiviral Res., 2007, 73, 228–231.
20. C. J. Woods, M. Malaisree, B. Long, S. McIntosh-Smith and A. J. Mulholland, 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.
21. F. S. Legge, A. Budi, H. Treutlein and I. Yarovsky, Protein flexibility: multiple molecular dynamics simulations of insulin chain B, Biophys. Chem., 2006, 119, 146–157.
22. J. P. Ma and M. Karplus, Molecular switch in signal transduction: reaction paths of the conformational changes in ras p21, Proc. Natl. Acad. Sci. U. S. A., 1997, 94, 11905–11910.
23. J. Apostolakis, P. Ferrara and A. Caflisch, Calculation of conformational transitions and barriers in solvated systems: application to the alanine dipeptide in water, J. Chem. Phys., 1999, 110, 2099–2108.
24. M. Karplus and A. Sali, Theoretical studies of protein-folding and unfolding, Curr. Opin. Struct. Biol., 1995, 5, 58–73.
25. L. Yang, G. Song, A. Carriquiry and R. L. Jernigan, Close correspondence between the motions from principal component analysis of multiple HIV-1 protease structures and elastic network modes, Structure, 2008, 16, 321–330.
26. A. Amadei, A. B. M. Linssen and H. J. C. Berendsen, Essential dynamics of proteins, Proteins: Struct., Funct., Genet., 1993, 17, 412–425.
27. D. M. F. vanAalten, J. B. C. Findlay, A. Amadei and H. J. C. Berendsen, Essential dynamics of the cellular retinol-binding protein – evidence for ligand-induced conformational changes, Protein Eng., 1995, 8, 1129–1135.
28. S. Bhakat, A. J. M. Martin and M. E. S. Soliman, An integrated molecular dynamics, principal component analysis and residue interaction network approach reveals the impact of M184V mutation on HIV reverse transcriptase resistance to lamivudine, Mol. BioSyst., 2014, 10, 2215–2228.
29. E. Laine, I. C. de Beauchene, D. Perahia, C. Auclair and L. Tchertanov, Mutation D816V Alters the Internal Structure and Dynamics of c-KIT Receptor Cytoplasmic Region: Implications for Dimerization and Activation Mechanisms, PLoS Comput. Biol., 2011, 7, 1–20.
30. E. F. Pettersen, T. D. Goddard, C. C. Huang, G. S. Couch, D. M. Greenblatt, E. C. Meng and T. E. Ferrin, UCSF chimera – a visualization system for exploratory research and analysis, J. Comput. Chem., 2004, 25, 1605–1612.
31. M. D. Hanwell, D. E. Curtis, D. C. Lonie, T. Vandermeersch, E. Zurek and G. R. Hutchison, Avogadro: an advanced semantic chemical editor, visualization, and analysis platform, J. Cheminf., 2012, 4, 17–33.
32. A. W. Goetz, M. J. Williamson, D. Xu, D. Poole, S. Le Grand and R. C. Walker, Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 1. Generalized Born, J. Chem. Theory Comput., 2012, 8, 1542–1555.
33. R. Salomon-Ferrer, A. W. Goetz, D. Poole, S. Le Grand and R. C. Walker, Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald, J. Chem. Theory Comput., 2013, 9, 3878–3888.
34. K. Lindorff-Larsen, S. Piana, K. Palmo, P. Maragakis, J. L. Klepeis, R. O. Dror and D. E. shaw, Improved side-chain torsion potentials for the Amber ff99SB protein force field, Proteins: Struct., Funct., Bioinf., 2010, 78, 1950–1958.
35. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys., 1983, 79, 926–935.
36. K. Kholmurodov, W. Smith, K. Yasuoka, T. Darden and T. Ebisuzaki, A smooth-particle mesh Ewald method for DL_POLY molecular dynamics simulation package on the Fujitsu VPP700, J. Comput. Chem., 2000, 21, 1187–1191.
37. P. Gonnet, P-SHAKE: a quadratically convergent SHAKE in O(n(2)), J. Comput. Phys., 2007, 220, 740–750.
38. S. Le Grand, A. W. Goetz and R. C. Walker, SPFP: Speed without compromise-a mixed precision model for GPU accelerated molecular dynamics simulations, Comput. Phys. Commun., 2013, 184, 374–380.
39. H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola and J. R. Haak, Molecular-dynamics with coupling to an external bath, J. Chem. Phys., 1984, 81, 3684–3690.
40. D. R. Roe and T. E. Cheatham III, PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data, J. Chem. Theory Comput., 2013, 9, 3084–3095.
41. P. A. Kollman, I. Massova, C. Reyes, B. Kuhn, S. H. Huo, L. Chong, M. Lee, T. Lee, Y. Duan, W. Wang, O. Donini, P. Cieplak, J. Srinivasan, D. A. Case and T. E. Cheatham, Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models, Acc. Chem. Res., 2000, 33, 889–897.
42. I. Massova and P. A. Kollman, Combined molecular mechanical and continuum solvent approach (MM-PBSA/GBSA) to predict ligand binding, Perspect. Drug Discovery Des., 2000, 18, 113–135.
43. V. Tsui and D. A. Case, Theory and applications of the generalized Born solvation model in macromolecular simulations, Biopolymers, 2000, 56, 275–291.
44. A. Onufriev, D. Bashford and D. A. Case, Modification of the generalized Born model suitable for macromolecules, J. Phys. Chem. B, 2000, 104, 3712–3720.
45. T. Hou, J. Wang, Y. Li and W. Wang, Assessing the Performance of the MM/PBSA and MM/GBSA Methods. 1. The Accuracy of Binding Free Energy Calculations Based on Molecular Dynamics Simulations, J. Chem. Inf. Model., 2011, 51, 69–82.
46. L. Xu, H. Sun, Y. Li, J. Wang and T. Hou, Assessing the Performance of MM/PBSA and MM/GBSA Methods. 3. The Impact of Force Fields and Ligand Charge Models, J. Phys. Chem. B, 2013, 117, 8408–8421.
47. D. R. Roe and T. E. Cheatham III, PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data, J. Chem. Theory Comput., 2013, 9, 3084–3095.
48. W. Humphrey, A. Dalke and K. Schulten, VMD: visual molecular dynamics, J. Mol. Graphics Modell., 1996, 14, 33–38.
49. A. Bakan, L. M. Meireles and I. Bahar, ProDy: Protein Dynamics Inferred from Theory and Experiments, Bioinformatics, 2011, 27, 1575–1577.
50. V. Karthick and K. Ramanathan, Insight into the Oseltamivir Resistance R292K Mutation in H5N1 Influenza Virus: A Molecular Docking and Molecular Dynamics Approach, Cell Biochem. Biophys., 2014, 68, 291–299.
51. A. C. Hurt, S. K. Leang, D. J. Speers, I. G. Barr and S. Maurer-Stroh, Mutations I117V and I117M and Oseltamivir Sensitivity of Pandemic (H1N1) 2009 Viruses, Emerging Infect. Dis., 2012, 18, 109–112.

---

Footnote

† Electronic supplementary information (ESI) available: RMSD vs. time and potential energy vs. time data for wild-type and I117V mutant of H5N1 neuraminidase complexed with oseltamivir from multiple MD trajectories is provided with supplementary materials available. See DOI: 10.1039/c4ra13494j

---