Introduction

The novel H7N9 influenza A virus (IAV) caused an outbreak of human infection in China in 2013. From March 2013 to May 2014, this IAV caused 440 human infections of which ~27 % were fatal [1]. Since then, this virus has additionally had 11 more confirmed human cases resulting in 4 more deaths on the 10 May 2016 (WHO). This new virus occurred by the multiple re-assortment of at least four avian influenza viruses; hemagglutinin subtype 7 (H7), neuraminidase subtype 9 (N9) and other six gene segments from two different groups of avian H9N2 viruses [2, 3]. The cause of poultry-to-human transmission is the mutation or genetic changes in domestic poultry, after which human can be exposed to infected poultry or contaminated environment in live poultry markets [2, 4–7]. Although there is no evidence for human-to-human transmission, the super re-assortment of this novel IVA is unpredictable in its emergence and incites the need for vaccine development [8, 9].

Many researches indicate that cross species transmission from poultry to humans is instigated by several mutations of proteins including the mutations in receptor binding site for hemagglutinin [10]. Such mutations include, the Q226L and G186V substitutions in hemagglutinin that result in stronger binding affinity to the human receptor [3], the mutation of a PB2 gene segment of E627K which is important for RNA replication in mammalian cell and the deletion of five amino acids (residue 69–73) in the stalk region of neuraminidase [2], which causes an enhancement in poultry to human adaptation [11].

Unlike in the case of the influenza virus where neuraminidase (NA) and the M2 ion-channel can be targeted, here with H7N9 influenza virus due to the S31N mutation only the neuraminidase can be targeted for effective treatment of the virus [2, 9, 12]. The four drugs that act as neuraminidase inhibitors (NAIs) are oseltamivir, zanamivir, peramivir and laninamivir (the latter of which is approved only in Japan). Neuraminidases are classified into two groups based on their primary sequence, except for N10 from bat-origin influenza-like virus and influenza B neuraminidase [13]. Group-1 NAs are N1, N4, N5 and N8 group-2 NAs are N2, N3, N6, N7 and N9 [14]. Whilst N1 and N2 are popular serotypes, generally found in epidemic and pandemic influenza viruses, this is a first time that N9 has infected human [15]. Sequence analysis of neuraminidase from some H7N9 strains have been shown to contain the R292K substitution after oseltamivir treatment such as in A/Shanghai/1/2013 [2, 16]. This R292K substitution is associated with multidrug resistance, particularly for the commonly used oseltamivir [15]. The arginine residue 292 (R292) is one of the arginine triad members (R118 and R371 in N2 numbering) located in the neuraminidase active site and directly interacts with the carboxylate group of sialic acid or with other drugs that bind in a similar fashion to sialic acid (see Fig. 1). In general, the R292K substitution, which causes the resistance of multi NAIs, is mostly found in the N2 and N9 neuraminidase subtypes [17, 18]. The novel H7N9 containing R292K mutation can be treated by NAIs, however, the protein- and virus-based assays clearly indicated drug resistance to oseltamivir and peramivir and a mild resistance to zanamivir and laninamivir, as listed in Table 1. Of particular importance is oseltamivir as it is used frequently in clinical treatment.

Fig. 1
figure 1

Structure of oseltamivir binding in the H7N9 neuraminidase active site a wild-type and b R292K substitution which results in multidrug resistance

Full size image
Table 1 Drug susceptibilities of plaque-purified H7N9 viruses and fold resistance of R292K substitution compare to those of wild-type
Full size table

Molecular dynamics (MD) simulation has been used successfully to elucidate the dynamics behaviors and equilibrium structures of biomolecules; whilst the calculated binding free energy predicts the strength of ligand binding affinity toward the protein target [19–22]. A combination of these two methods could result in better understanding of drug binding mechanism and source of drug resistance, which the obtained information is further useful for designing and developing the more potent inhibitors [23–25]. Although MD simulations of drug-H7N9 neuraminidase complexes have previously been carried out by others and their binding affinities predicted [26, 27], their initial structures were prepared based on homology modeling. Since then the crystal structure of H7N9 neuraminidases have been made available. This work aims to understand oseltamivir binding to the novel H7N9 neuraminidase and provide insight into the role of the R292K mutation at the molecular level by using MD simulations and binding free energy calculations of molecular mechanics combined with Poisson–Boltzmann (PB) or generalized Born (GB) electrostatics solvation and solvent accessible surface area nonpolar solvation (SASA); MM/PB(GB)SA methods. These well-known binding free energy approaches have been successful applied to various ligand binding affinities and drug resistances [28–42].

Materials and methods

Complex preparation

Oseltamivir binds to both wild-type and R292K mutation H7N9 neuraminidase. Complexes were prepared from crystal structure of the unbound neuraminidase Protein Data Bank (PDB) structures, with codes of 4MWJ and 4MWL respectively, with superimposition to their oseltamivir structure of 4MWQ and 4MWW respectively [15]. These systems were then prepared for calculation and carried out using the AMBER 14 software package [46]. The partial charges and parameters of oseltamivir was obtained from previous work [47], and the force fields of the protein were set to the ff12SB model of Amber 14. The ionized states of all side chains of lysine (K), arginine (R), aspartate (D), glutamate (E) and histidine (H) were configured to be at pH 7.0 using PROPKA3.1 [48]. All disulfide bonds between pairs of cysteine (C) residues were defined to stabilize the protein structure. The LEaP module added the missing hydrogen atoms in the oseltamivir-NA complex. A crystal calcium ion was kept, while water molecules were deleted. The placevent algorithm based on the 3D-RISM distribution function was employed to firstly solvate TIP3P water model around solute molecule [49]. The TIP3P waters were then filled with a minimum distance of 10 Å from the solute to the edge of the simulation box (88 × 90 × 81 Å3) and two Cl− ions were added to neutralize the system.

MD simulation and protocol

To optimize the initial structure before performing MD simulations, the complexes were minimized in the following order, hydrogen atoms, water molecules and then the entire system implementing a 1000 steps using the steepest descent algorithm followed by 2000 steps with the conjugated gradient algorithm. The energy minimization process and MD simulations were performed using the pmemd. CUDA module of AMBER. All bonds involving hydrogen atoms were constrained using the SHAKE algorithm [50]. To reduce the computational time of non-bonded interactions, a cutoff distance of 10 Å was used and the particle mesh Ewald method [51] was employed to calculate the long-range electrostatic interaction.

The simulation employed the periodic boundary condition using Langevin dynamics with a collision frequency of 2 ps−1 and a simulation time step of 2 fs. The system was heated up from 10 to 300 K for 100 ps using canonical ensemble (NVT) in the thermalization step. Subsequently, the equilibration step was firstly conducted for 100 ps of NVT followed by NPT ensemble at 1 atm and 300 K until 30 ns. Finally, the system was further simulated for 50 ns. The production phase of the root mean square displacement (RMSD) analysis was assigned from 30 to 50 ns.

MM/PB(GB)SA binding free energy calculations

The binding free energy is the free energy difference of three reactants; protein, ligand and complex, according to Eq. 1. The binding free energy between protein and ligand is the summation of free energy change during ligand binding \((\Delta G_{MM}^{{}} )\) and the solvation free energy \((\Delta \Delta G_{solv}^{{}} )\).

$$\Delta G_{bind} = G_{complex} - G_{receptor} - G_{ligand} = \Delta G_{MM} + \Delta \Delta G_{solv}$$
(1)
$$\Delta G_{MM}^{{}} = E_{MM} - T\Delta S = (\Delta E_{MM}^{ele} + \Delta E_{MM}^{vdw} ) - T\Delta S$$
(2)
$$\Delta \Delta G_{solv}^{{}} = \Delta \Delta G^{polar}_{solv} + \Delta \Delta G^{nonpolar}_{solv}$$
(3)
$$\Delta \Delta G_{solv}^{nonpolar} = \gamma \Delta SASA$$
(4)

The molecular mechanical energy is the combination of electrostatic and van der Waals interactions as outlined in Eq. 2, similar to the solvation free energy it can be decomposed into electrostatic and nonpolar contributions as showed in Eq. 3. The MM/PBSA and MM/GBSA methods [32, 52–60] compute the polar solvation free energy by Poisson–Boltzmann or generalized Born models, respectively, whereas their nonpolar solvation energy was obtained by solvent accessible surface area (SASA) using a probe radius of 1.4 Å and a surface tension constant (γ) of 0.0072 kcal mol−1 Å2 according to Eq. 4. The binding free energy calculation was averaged over 100 snapshots extracted from the production phase using the MMPBSA.py program [61]. The electrostatic solvation free energies are calculated with a manner of each method. MM/PBSA computes the polar contribution by solving the linear Poisson–Boltzmann equation using the pbsa module implemented in sander. The grid spacing was set to 0.125 Å and a probe radius of water molecule is of 1.4 Å. In contrast for MM/GBSA, the polar solvation free energy is carried out with the Generalized Born (GB) model developed by Onufriev and colleagues (igb = 2) [62] used in this study. The dielectric constants for solute and solvent of 1 and 80 are applied. To reduce the high computer consuming, the structural entropy calculation is estimated over the 25 snapshots of the MD structures by normal mode analysis in the gas phase [63]. Each structure was entirely minimized for 1000 steps until reaching a local minimum by considering energy gradient with less than 0.01 kcal mol−1 Å−1. The Hessian matrix is diagonalized to obtain the vibrational frequencies, while the translational and vibrational entropies are computed by statistical mechanics formula [53, 64].

Results and discussion

  • Stability of complexes

The stabilities of each complex were evaluated using root-mean square displacement (RMSD) to explore the maintained fluctuation of the MD structure compared to their initial structure. In Fig. 2, the RMSD was calculated in terms of backbone atoms, which belong to the neuraminidase protein and the ligand atoms of oseltamivir versus simulation time. For backbone RMSD, both wild-type and R292K complexes show increasing RMSD from ~0.5 to 1.0 Å during the first 10 ns and after which it oscillates at this value until 30 ns. Finally, the fluctuation stabilizes after 30 ns. With regards to the ligand RMSD of oseltamivir the values varied more ranging from 0.5 to 2.0 Å for the entire simulation. It can be clearly seen that oseltamivir structure of R292K system has significantly higher fluctuations compared to the wild-type system. As a result, the MD structures from 30 to 50 ns were adopted as the production phase for further analysis.

Fig. 2
figure 2

RMSD calculations in term of the backbone atoms of the neuraminidase protein (upper) and the oseltamivir ligand atoms (lower) versus the simulation time

Full size image

  • Binding free energy calculation of drug susceptibility

The aspects of this study are to understand the oseltamivir efficiency towards the H7N9 neuraminidases and to explore the source of reduce susceptibility concerning the R292K mutation. The MM/PBSA and MM/GBSA binding free energies and their components are shown in Table 2 and are also compared to the experimental values.

Table 2 Binding free energy of oseltamivir toward wild-type and R292K neuraminidases using MM/PB(GB)SA methods
Full size table

In the wild-type complex, MM/PBSA method qualitatively predicts the binding free energy of −14.7 kcal mol−1, which is consistent with the experimental values of −12 to −13 kcal mol−1 [15, 43–45], while MM/GBSA approach underestimates the value with a prediction of −28.5 kcal mol−1. The total binding free energy (\(\Delta G_{bind}\)) of these methods is a subtle balance between the negative value of the binding interaction (\(\Delta G_{MM}\)) and the positive value of the solvation free energy (\(\Delta \Delta G_{solv}\)) called the dehydration penalty which involves the loss of solvation around the protein after ligand binding and is clarified in more detail in a previous study [65]. The electrostatic contribution plays a major role for oseltamivir binding as can be seen from the large free energy values of the electrostatic interaction \(\Delta E_{ele}^{MM}\) and the polar solvation free energy \(\Delta \Delta G_{solv}^{polar}\). This is due to the zwitter ion nature of oseltamivir accommodating the highly charged binding pocket of neuraminidase.

The increase in magnitude of binding free energy in the R292K mutated system was calculated for MM/PBSA (−4.5 kcal mol−1) and MM/GBSA (−19.8 kcal mol−1) to be 10.2 and 8.7 kcal mol−1, respectively, compared to those of wild-type. Similar to the wild-type complex, the MM/GBSA method underestimates the prediction of the binding free energy compared to the MM/PBSA approach (experimental values of −5.2 or −7.2 kcal mol-1). The source of reduced oseltamivir susceptibility conferred by the R292K substitution is a reduction in binding interaction (negative value) and solvation free energy (positive energy). Both methodologies show a loss of drug interaction by 18.5 kcal mol−1, and a decrease in dehydration penalty resulting in a change of solvation energy by 8.3 kcal mol−1 for MM/PBSA and 9.8 kcal mol−1 for MM/GBSA. In the following sections we will further elucidate the effect of the R292K substitution towards oseltamivir binding.

  • Change of drug binding interaction

The favorable electrostatic interaction of oseltamivir binding is obtained from strong hydrogen bond (H-bond) formations between polar moieties of the drug (including –COO−, –NHAc and –NH3 +) and the highly charged binding pocket of neuraminidase. In addition, the nonpolar side chain of the pentyl group gives van der Waals contribution stabilizing the hydrophobic pocket formed by the E276 and R224 residues. The H-bond formation is illustrated as a percentage occupation, where occupation is defined as having maximum distance of 3.5 Å between hydrogen donor (D) and acceptor (A) together with the minimum bond angle of 120° between three atoms (D…H…A) depicted in Fig. 3 for the wild-type and the R292K mutant neuraminidases.

Fig. 3
figure 3

Percentage occupation of the hydrogen bond interactions between oseltamivir function groups and the H7N9 neuraminidase binding residues for a wild-type and b R292K strains. Also shown is a close up of oseltamivir in the neuraminidase binding pocket of c wild-type and d R292K strain taken from the last snapshot

Full size image

In Fig. 3, oseltamivir forms H-bond interaction with residues of the arginine triad (R118, R292 and R371, according to N2 numbering), R152, E119 and D151 as is generally found to be the case in other neuraminidase subtypes [21, 33, 47, 65]. In the R292K complex, there is a complete loss of the carboxyalte-K292 interaction as previously found in the H1N1 subtype [21], however, this is compensated for by the R118 interaction (Fig. 3b). A smaller side chain of lysine (K) compared to arginine (R) at position 292 cannot form the necessary H-bond with the carboxylate group (distance of around 5 Å), leading to the larger cavity and shifting toward the R118 residue of ligand (Fig. 3d). It is worth noting that the conformation of this K292 mutated residue does not perturb the formation of hydrophobic pocket made by E276 and R244 residues as found in the crystal structure [15, 27].

The conformational change of the pentyl group was evaluated by the distribution of torsion angles (Fig. 4). Both wild-type and R292K systems have different angular distributions, particularly for τz. The torsion distribution of τz in the mutant system dramatically changes from the normal ~60° of the wild-type, while τx and τy show the same pattern but different distributions. This suggests that the drug in the R292K complex is oscillating. This oscillation can clearly be seen from the high fluctuation of ligand RMSD during simulation (in Fig. 2). The twist of pentyl group was found to disrupt the hydrophobic binding in a previous study [26].

Fig. 4
figure 4

Distribution of the torsion angles of the oseltamivir pentyl group calculated while in complexation with the binding pockets

Full size image

To summarize, the unsuitable accommodation of oseltamivir in R292K neuraminidase results in a decrease of drug binding interaction by 13.8 kcal mol−1 for electrostatic interaction and 5.1 kcal mol−1 for van der Waals interaction compared to the wild-type complex.

  • Change of solvation conferring R292K substitution

The latter source of the high-level oseltamivir resistance toward R292K substitution is due to the decrease of the positive contribution of the solvation free energy resulting in a smaller dehydration penalty compared to the wild-type. As mentioned in previous sections, the shorter side chain of K292 causes poor binding with the drug carboxylate group and a larger cavity. This results in greater accessibility for water molecules towards the neuraminidase binding pocket as evident by the radial distribution function (RDF) around oseltamivir heteroatoms and residue 292 (lables shown in Fig. 3c, d) depicted in Fig. 5.

Fig. 5
figure 5

Radial distribution function g(r) of water atom and the integrated up coordination number around a all heteroatoms of oseltamivir and b the residue of 292

Full size image

In Fig. 5a, the first peak of water distribution around oseltamivir heteroatoms occurs at 3.5 Å and is mostly sharp, which suggests that the accessible water molecules strongly interact with oseltamivir. The nonzero value at the terminal of the first peak as found in the O1, O2 and N1 graphs indicate a high degree of transfer of water molecules in this shell. Switching of water molecules around the oseltamivir molecule was found to occur between the O1 and O2 atoms at the oseltamivir carboxylate group (see illustration in Fig. 3c, d), while the rest of the heteroatoms are conserved for both wild-type and R292K systems. Additionally two water molecules were found to bind around the K292 residue, compared to the one molecule of the wild-type (Fig. 5b). This mutated residue was strongly solvated by three water molecules as can be seen by the approximately close to zero values at the fist minimum, which is in good agreement with the crystal structures of H7N9 and H11N9 [15, 66]. Thus, the R292K mutation of neuraminidase can affect the accessibility of water in the binding pocket, which causes a decreased dehydration penalty.

Conclusion

In this study, oseltamivir binding towards the novel H7N9 influenza neuraminidase and the high resistance R292K substitution was investigated using MD simulations and MM/PB(GB)SA binding free energy calculations. MM/PB(GB)SA methods were able to predict the susceptibility of oseltamivir toward both the wild-type and the R292K strain and is in good agreement with the energy converted from the IC 50 data, of which the MM/PBSA approach was shown to give a more reliable prediction. A structural entropy calculation was necessary to improve the qualitative of the predicted binding free energy. Similar to the previous study of influenza B neuraminidase, the main factor for oseltamivir binding toward H7N9 neuraminidase is a subtle balance of the large negative value of binding interaction and the large positive contribution of solvation conferring the dehydration penalty during drug binding. The major interaction of oseltamivir binding is the electrostatic contribution, as is generally found for all neuraminidase subtypes, caused by the opposite charges of oseltamivir and neuraminidase binding residues. A dramatic increase in magnitude of binding free energy was found in the R292K system compared to the wild-type. The role of R292K mutation in H7N9 neuraminidase is to decrease the binding interaction as was shown by the high fluctuation of drug RMSD in the mutant pocket and the reduction of the dehydration penalty resulting in an increased accessibility for water molecules around the K292 mutated residue and the carboxylate group of drug. The results from this work suggest that a hydrophilic bulky group as found in zanamivir and laninamivir, instead of the pentyl group of oseltamivir, promises to be a potent neuraminidase inhibitor by maintaining a strong interaction toward residue 292 of both arginine and lysine amino acids.