Computational study of binding of oseltamivir to neuraminidase mutants of influenza A virus

Oseltamivir (OTV), which targets the neuraminidase (NA) of Influenza A virus (IAV), has been reported to develop resistance. Here, we performed a computational study on the binding modes of OTV in the wild-type and popular mutants of IAV NA (E119A, E119D, E119G, H274Y, I117T, I117V, I117V-E119A, K150N, N294S, R292K, V116A


INTRODUCTION
Influenza, caused by influenza viruses, is an acute respiratory infection of a potentially deadly disease that infects hundreds of millions of people each year around the globe [1].The burden of influenza exists despite available vaccines and with the current coronavirus disease 2019 infections, which complicates health care systems and indicates the urgent need for effective anti-influenza drugs.
Influenza virus consists of four genera, i.e., types A, B, C, and D [2], in which the Influenza A virus (IAV) is the most pathogenic to humans among the four genera.Oseltamivir (OTV), also known as Tamiflu, is an antiviral medication that is commonly used to treat and prevent influenza A and B infections.It blocks the activity of NA, which allows the virus to spread from infected to healthy cells.It can help reduce the severity and duration of the illness, as well as the risk of complications such as pneumonia.
Mihajlovic and Mitrasinovic [11] studied mutations at His274 which slightly enhanced or reduced the susceptibility of OTV to NA.The effect of H274Y mutation in drug-target interactions was also studied by Karthick and Ramanathan
Binding affinities calculated were −12.857, −7.760, −11.683, −7.860, −10.575, −10.887, −10.898, −7.647, −10.999, −10.947, −10.534, −10.825, and −11.080 kcal/mol, for the WT, E119A, E119D, E119G, H274Y, I117T, I117V, I117V-E119A, K150N, N294S, R292K, V116A, Y252H systems, respectively (Table S1).The mutation reduced the docked binding energy of [12], which implied the stable binding of the OTV with wildtype (WT) NA as compared to mutant-type NA.In both studies, molecular docking and short MDS were used.Nguyen et al. [13] used the smooth reaction path generation method to calculate binding free energies and found that the binding free energies of NA A/H5N1 in WT and H274Y and N294S mutants with OTV show good agreement with experimental results.To the best of our knowledge, there is no previous study on OTV resistance that covers more comprehensive residue mutations and longer MDSs than our present work.The present study would reveal the atomistic details about designing OTV analogs having better activity toward the NA of IAV.

COMPUTATIONAL METHOD
Homology modeling, molecular docking (Glide XP), MDS, EDA, and molecular mechanics with generalized born and surface area solvation (MM-GBSA) calculations were used to study the binding of OTV to WT and mutant proteins as described in other papers [14,15].

MDS, analysis, MM-GBSA binding energy calculations
Each protein-ligand system was used to construct MDS systems using the OPLS_2005 force field [19].Na+ ions were added to a NaCl concentration of 0.15 M NaCl to neutralize the systems and were used to model the protein-ligand systems.All MDSs were performed for 200 ns using the Desmond simulation package [20].
The system relaxation was carried out based on previous work [17,[21][22][23].The analysis of simulation results was computed over the entire trajectory using the Desmond trajectory clustering tool [17,24].
MM-GBSA binding energies were predicted for each protein-ligand system as described in the previous works [10,17,25].

Energy decomposition analyses
The EDA calculations were conducted using the absolutely localized molecular orbitals (ALMO)-EDA scheme with the ωB97M-V functional and the Def2-SVPD basis set [26][27][28] using the Q-Chem program suite [29].This analysis was performed for mutations and residues for which the proteinligand interaction analysis resulted in interactions above 50% (vide supra).The EDA analyses allowed us to measure the contributions of dispersion (∆E DISP ), polarization (∆E POL ), and charge transfer (∆E CT ) interactions to the overall intermolecular interaction energies.OTV to NA.The 2D interactions of OTV and WT and mutant systems are shown in Figure S1.higher than that in I117T (1.66 Å), H274Y (1.74 Å), and N294S (1.78 Å).

Cluster analysis
H-bonds with Glu119 and Asp151 were observed in almost all mutant systems except for E119A, E119G, E117V-E119A, R292K, and V116A, whereas H-bonds with Glu277 H294S, and V116A systems.Other interactions with Arg118 and Arg371 were noted in the K150N and V116A systems, respectively.Figure 4 shows the 2D and 3D interactions of WT OTV taken from the most populated cluster, while those of mutants systems were described in Figure S2.

EDA results
The EDA calculation revealed the nature of the H-bonding and dispersion interactions between OTV and the active site residues.Each of the small models comprises a drug in addition to a key residue of the active site.The geometries of the dimers were obtained from converged MD trajectories.Table 3      interaction energies between OTV and key residues in the active site.

Continued
Inspection of Table 3 reveals that, with a few exceptions, the polarization interaction is the dominant interaction, followed by charge-transfer and dispersion interactions.It is instructive to begin with a discussion on the WT system.The strongest overall interaction was obtained for the ARG371 residue, where the sum of the ∆E DISP , ∆E POL , and ∆E CT components (∆E TOT ) amounted to as much as −188.3kJ/ mol.We note that polarization interactions account for −99.5 kJ/ mol of this total interaction.The ARG292 residue also interacts strongly with the drug, with an overall interaction energy of −161.6 kJ/mol.For both ARG371 and ARG292 residues, we obtained high percentage interactions of 95%-98% from the protein-ligand interaction analysis.We also obtain a strong overall interaction of −123.9 kJ/mol with the ASP151 residue, where the smaller overall interaction energy is mainly due to a reduction in the polarization component to merely −55.8 kJ/mol (Table 1).Accordingly, the protein-ligand interaction analysis resulted in a smaller value of 87%.Finally, for the ARG118 residue, we obtained a low overall interaction of −17.0 kJ/mol, consistent with an even smaller percentage interaction of 62% obtained from the protein-ligand interaction analysis.
Let us now discuss the E119A mutant.In this mutation, we observed a strong interaction with the ARG292 residue, where the overall interaction amounts to −157.9 kJ/ mol.Polarization interaction accounts for approximately half of the total interaction (kJ/mol).We also obtained a strong overall interaction with the ASP151 residue at −120.7 kJ/mol.However, dispersion interactions play a lesser role.These results are consistent with the protein-ligand interaction analysis, where both residues resulted in percentage interactions between 81% and 91%.For the related E119D mutation, we obtained a significant overall interaction of −137.0 kJ/mol with the ASP151 residue.Again, the polarization component accounted for nearly 50% of the overall interaction energy.It should be noted that the protein-ligand interaction analysis resulted in 99% interaction for this residue.
For the H274Y mutation, we obtained a very large overall interaction of −189.5 kJ/mol with the ARG371 residue, where the polarization component accounts for over 50% of the overall interaction energy, namely −99.1 kJ/mol.Protein-ligand interaction analysis resulted in 96% interaction for this residue.For the GLU119 residue, we obtain a smaller overall interaction of −85.9 kJ/mol, consistent with a percent interaction of 86%.For the ASP151 residue, we obtained a low overall interaction of merely −12.8 kJ/mol, consistent with an interaction of only 50%.
For the I117T mutation, we obtained a significant interaction of −124.5 kJ/mol with the GLU119 residue.The interaction with the ASP151 residue was significantly smaller and amounted to −52.8 kJ/mol.Both interactions are dominated by the polarization component, and both are consistent with the percent interactions of 86% and 79% obtained from the protein-ligand interaction analysis.For the I117V mutation, we obtained similar interaction energies ranging between −97.2 and −111.4 for the three residues GLU119, ASP151, and ARG292.
For the K150N mutation, the interaction energies of the three residues varied widely.For example, we obtain the following interaction energies −133.8 (ARG118), −77.8 (ARG292), and −14.8 (ASP151).Therefore, ARG118 plays a significantly more important role in drug binding.We note, however, that in contrast to the other mutations, for the K150N mutation, there was no agreement with the proteinligand interaction analysis.For example, we obtained percent interactions of 51% (ARG118), 54% (ARG292), and 87% (ASP151).
In the case of the N294S mutation, three residues interacted strongly with the drug.ARG292, ARG371, and GLU119 resulted in binding energies ranging between −124.0 and −129.0 kJ/mol.For the R292K mutation, we obtain a similarly strong interaction of −122.2 kJ/mol with the GLU119 residue.For the V116A mutation, both arginine residues ARG292 and ARG371 interact strongly with the drug, with interaction energies of −152.5 and −121.4 kJ/mol, respectively.Finally, for the Y252H mutation, the ASP151 and GLU119 residues result in appreciable interaction energies of −104.1 and −119.1 kJ/mol, respectively.

The root mean square fluctuation (RMSF) values
Figure 5 depicts the RMSF of the protein Cα atoms in the WT and mutant strains.The highest peaks were for residues of the carboxyl and amino ends of the protein.The high peaks around 29, 60, 167, 188, 329, and 351 were protein loops, while the key residues interacting with the protein were observed to be stable.R292K (1.57Å), V116A (0.98 Å), and Y252H (1.382 Å), but it was higher than those in E119D (0.76 Å) systems, which indicated that OTV was more stable in the WT system than in I117V-E119A but it was less stable compared to in E119D systems.The atom numbers 16-20 (ethylpropoxy group) in I117V-E119A were considered the most fluctuant atoms.

Ligand RMSF values
The OTV dihedral angle profiles Figure 7 depicts the dihedral angle of OTV, in which the rotatable bonds colored light blue, light green, yellow, light purple, red, and orange show a wider distribution in mutant systems.Wider dihedral angle distributions were also observed in yellow for all mutants and light green in the E119G, I117T, I117V-E119A, and V116A systems.Briefly, the mutant systems induced a wider distribution of dihedral angles within OTV.
OTV is an FDA-approved drug that inhibits the NA activity of IVA.Previous experimental studies have shown that Arg118, Arg292, and Arg371 play vital roles in the active site of NA, whereas Glu276 interacts with acetamido from the substrate to form hydrogen bonds.Asp151 and Glu119 are responsible for cavity widening at the active site of NA [30].These experimental observations correlated well with the results of the present study for the WT system.
Recent findings suggest that resistance to OTV arises from a NA mutation in IAV.In vitro assays revealed mutations, including E119A, E119D, E119G, H274Y, I117T, I117V, I117V-E119A, K150N, N294S, R292K, V116A, and Y252H, resulting in decreased sensitivity to OTV.Here, we show that NA mutations have reduced OTV affinity due to reduced electrostatic interactions that become more positive in all systems.In the K150N system, the mutation induced a more positive van der Waals interaction, resulting in weakened binding of OTV to IAV.Meanwhile, the I117V-E119A and R292K systems also exhibited low OTV binding energies owing to increased electrostatic values during the simulation.In addition, it may be triggered by reduced occupancy and loss of interaction with several key residues such as Arg292, Arg371, Glu119, and Asp151.

CONCLUSION
The molecular details of OTV binding to IAV NA in WT and mutant systems were explored using a computational study.The key residues that interact with OTV include Arg118, Arg292, Arg371, Arg152, Glu276, Asp151, and Glu119, which are in line with the experimental results.The binding modes of OTV changed, resulting in reduced affinity in all

Figure 1 .
Figure 1.The NA is complex with OTV (orange color).The mutated residues are shown.

Figure 4 .
Figure 4.The active site conformation and 2D interactions of the WT-OTV system taken from the most populated cluster.H-bonds were formed between the drug and Glu119, Asp151, Arg292, and Arg371.

Figure 6 .
Figure 6.The RMSF values for atoms of OTV.

Figure 7 .
Figure 7.The dihedral angle of OTV profiles during 200 ns MDS which is the conformational progression of the nine rotatable bonds of OTV.The dial plots describe the conformation of the torsion throughout the course of the simulation.The beginning of the simulation is in the center of the radial plot and the time evolution is plotted radially outward.The bar plots summarize the data on the dial plots, by showing the probability density of the torsion.In addition, the 2D structure of OTV in the top panel is for reference.

Figure 7 .Figure 7 .
Figure 7.The dihedral angle of OTV profiles during 200 ns MDS which is the conformational progression of the nine rotatable bonds of OTV.The dial plots describe the conformation of the torsion throughout the course of the simulation.The beginning of the simulation is in the center of the radial plot and the time evolution is plotted radially outward.The bar plots summarize the data on the dial plots, by showing the probability density of the torsion.In addition, the 2D structure of OTV in the top panel is for reference.Continued

Figure 7 .
Figure 7.The dihedral angle of OTV profiles during 200 ns MDS which is the conformational progression of the nine rotatable bonds of OTV.The dial plots describe the conformation of the torsion throughout the course of the simulation.The beginning of the simulation is in the center of the radial plot and the time evolution is plotted radially outward.The bar plots summarize the data on the dial plots, by showing the probability density of the torsion.In addition, the 2D structure of OTV in the top panel is for reference.

Table 2 .
The binding energies calculated for the last 100 ns of each trajectory.

Table 3 .
Breakdowns of the dispersion (∆E DISP ), polarization (∆E POL ), and charge transfer (∆E CT ) interaction energies between OTV and key residues in the active site obtained from second-generation ALMO-EDA at the _B97M-V/Def2-SVPD level of theory (in kJ/mol).