In silico modification of oseltamivir as neuraminidase inhibitor of influenza A virus subtype H1N1.

This research focused on the modification of the functional groups of oseltamivir as neuraminidase inhibitor against influenza A virus subtype H1N1. Interactions of three of the best ligands were evaluated in the hydrated state using molecular dynamics simulation at two different temperatures. The docking result showed that AD3BF2D ligand (N-[(1S,6R)-5-amino-5-{[(2R,3S,4S)-3,4-dihydroxy-4-(hydroxymethyl) tetrahydrofuran-2-yl]oxy}-4-formylcyclohex-3-en-1-yl]acetamide-3-(1-ethylpropoxy)-1-cyclohexene-1-carboxylate) had better binding energy values than standard oseltamivir. AD3BF2D had several interactions, including hydrogen bonds, with the residues in the catalytic site of neuraminidase as identified by molecular dynamics simulation. The results showed that AD3BF2D ligand can be used as a good candidate for neuraminidase inhibitor to cope with influenza A virus subtype H1N1.

Introduction

Neuraminidase has become the main target for drug design against influenza virus due to its highly conserved catalytic site and its essential role in influenza virus replication[1]-[5]. Oseltamivir is an antiviral drug against neuraminidase that is useful for preventing viral replication in the last stage of the viral live cycle[6]. It directly interacts with the catalytic residues of the neuraminidase catalytic site, while the framework residues stabilize the enzyme structure[7]. Resistance to oseltamivir has already become a widespread phenomenon, and prediction of the neuraminidase structure shows that the resistance could be several times higher than resistance to zanamivir, a neuraminidase inhibitor[8],[9]. Mutations at the conserved residues of neuraminidase appear to associate with oseltamivir resistance in a subtype specific manner[6]. Arg292Lys, Asn294Ser, and His274Tyr are three known types of mutations that manifest oseltamivir resistance in H1N1 influenza virus[1],[9]-[11]. The frequency of His274Tyr mutant viruses is 64%, although it varies among different countries[9],[12] . A high prevalence of the mutant virus is seen in South Africa, New Caledonia, New Zealand, the Philippines, and Australia[13], modest in Singapore, Malaysia and Thailand, but insignificant in Macau and Taiwan[12].

The search for effective drugs that can cope with H1N1 is still ongoing. Li et al. searched zinc fragment database and found the lead compound of Neo6 as a feasible drug candidate[14]. Drug resistance tendecy of oseltamivir towards H5N1 provided valuable data for H1N1 drug design[1],[11]. Six analog inhibitors as drug candidates against H5N1 were suggested by Du et al.[15]. A more robust binding energy than that of zanamivir and oseltamivir was found in compound AG7088 and its derivatives[16],[17].

Oseltamivir has good bioavaibility and is effective for treatment and prevention of epidemic influenza infection in adults, adolescents and children (≥ one year of age)[18]. Modifcation of oseltamivir to cope with influenza A virus subtype H1N1 has been attempted[1],[11], which is based on the properties of amino acid residues in the catalytic site of neuraminidase[19]. Drug-target interactions are being investigated using structural bioinformatics tools[20],[21]. In this study, we focused on the modification of the oseltamivir functional groups and molecular docking was performed to determine the best candidate ligand based on the lowest binding energy and interaction. Molecular dynamics simulation is a useful theoretical tool for analyzing the protein-ligand interactions with atomic resolutions based on classical mechanics[22]-[25]. We also carried out molecular dynamics simulation to evaluate the interaction of candidate ligands and the enzyme in the hydrated state at two different temperatures.

Materials and methods

Sequence alignment and homology modeling of neuraminidase

Neuraminidase sequences with mutations were downloaded from National Center for Biotechnology Information (NCBI) influenza virus sequence database (http://www.ncbi.mlm.nih.gov/genomes/flu/). The multiple sequence alignment method was based on the CLUSTAL W2 program (www.ebi.ac.uk/Tools/clustalw2/index.html)[11],[26]. The alignment result was forwarded into homology pipeline for further processing of 3D structure determination. The homology modeling was performed using the Swiss model, which can be accessed through http://www.swissmodel.expasy.org/SWISS-MODEL.html. 3CKZ chain A [Protein Data Bank (PDB) code] with mutation His274Tyr as template was applied to build the latest structure of N1. Validation of the 3D structure from the homology modeling was performed by protein geometry program and superimposed by superpose program in MOE 2009.10 software[9]. Based on the superimposition, the root-mean-square deviation (RMSD) was calculated to identify structural similarity between template model mutated with 3D structure from homology modeling.

Ligand preparation of oseltamivir modification

Modification of the oseltamivir functional groups were conducted with alcohol, aldehyde, ketone, carboxylic acid, esther, amide, O-glycoside[27] and apioside group[28] using ACD Labs software[29]. Modified oseltamivir was built into 3D structure using ACD Labs software. 3D shape was obtained by saving it in the 3D viewer of ACD Labs. Furthermore, the output format was converted into Molfile MDL Mol format using software Vegazz to conform the docking process[30]. Ligand was washed with a computer program; adjustments were made with the ligand partial charge and partial charge optimization using MMFF94x forcefield[9],[31]. The conformation structure energy of ligands was minimized using the root-mean-square (RMS) gradient energy with 0.001 kcal/mol[9]. Other parameters were in accordance with existing default in the MOE 2009.10 software[32].

Molecular docking

MOE 2009.10 was applied for optimization and minimization of the 3D structure of the enzyme with the addition of hydrogen atoms[9]. Protonation was employed with Protonate 3D programs to introduce changes in ionization step. The ‘wash’ process was started by optimization pipeline, on MOE database viewer. Furthermore, partial charges and force field modification were applied with MMFF94x[9],[30]. Solvation of enzymes was performed in the form of a gas phase with a fixed charge, RMS gradient of 0.05 kcal/ mol, and other parameters by using the standard in MOE 2009.10 software[9],[30],[33].

The initialization of the docking process was applied with MOE 2009.10[9]. Docking simulations were performed by Compute-Simulation dock program. During the docking procedure, ligands were flexible, whereas the receptor was fixed. The application of triangle matcher for the placement method with the iteration of total 1,000,000 energy readings per position and parameter tuning was performed with the standard default value of MOE[9],[33]. A total of 1,000 populations for scoring function were utilized for LondonDG and repetition force field[9],[30],[33]. The first 100 repetitions and the second setting were shown only as one of the best results.

Scoring using London DG estimated Gibbs binding energy (ΔGbinding) from pose of ligand to enzyme was based on the calculation:Where, c  =  rotation entropy rate and translation that can be obtained or released, E  =  flexibility of ligand energy, f  =  imperfection size geometry of hydrogen bonds, c  =  energy of an ideal hydrogen bond, f  =  imperfection size geometry of metal ligations, cM  =  ideal metal ligation energy and D  =  atom desolvation energy[31],[34]-[36]. The analysis of the ligand interactions with enzyme at both temperatures of 300 K and 312 K was employed by LigX tools in the MOE 2009.10 software.

Molecular dynamics simulation

Molecular dynamics simulation was carried out at 2 different temperatures using Simulations-Dynamics program in MOE 2009.10 software[37]. The solvation in potential setup was born with RMS gradient of 0.05[38]. Geometry optimization and energy minimization of 3D structure neuraminidase complex ligands were perfomed with MOE 2009.10 software[31]. Geometry optimization using partial charge was current force fields. Energy minimization was employed with MMFF94x. Ensemble parameter was performed: NVT (N, number of atoms; V, volume; T, temperature), NPA algorithm, and Cutoff restraint 6Å. Dynamics simulation stage involved initialization for 40 ps, main simulations as equilibrium and production stage for produced trajectory.

Molecular dynamics simulation of complex enzyme-ligand at temperature of 300 K was employed with main simulation for 5,000 ps and cooling for 20 ps[39]. Molecular dynamics simulation at temperature of 312 K was employed from heating at temperature of 300 K to 312 K for 20 ps. Main simulation was employed for 5000 ps and cooling for 20 ps. Position, velocity and acceleration were saved every 0.5 ps. The other parameters were performed in accordance with default MOE dynamics parameters.

Born solvation caused the solvation energy (Esol) to be calculated in potential energy molecular system function from atom coordinate:Where, Estr  =  stretching energy, Eang  =  angular energy, Estb  =  stretching and bending energy, Eoop  =  operating energy, Etor  =  energy torsional, Evdw  =  intermolecules van der Waals energy, Eele  =  electrostatic energy, Esol  =  solvation energy, dan E  =  residual moving energy so that we can see clearly Esol in the equation[39]–[41].

Interaction between ligand and residues after simulation was analyzed using LigX Ligand Interaction in MOE 2009.10 software[42]. Visualization of different conformations that occured during simulation was analyzed using Surface and Map program in MOE 2009.10 software.

Results

Homology modeling

Neuraminidase sequence > gi|237651250|gb|ACR08499.1| neuraminidase [Influenza A virus (A/Auckland/1/2009(H1N1))] was used as the target sequence and had 91.123% similarity with the template of 3CKZ chain A (PDB code). Accordingly, the alignment results were employed to build homology model for N1 and then superimposed with the template[1]. Homology modeling for neuraminidase of influenza A virus H1N1 was obtained using template 3CKZ chain A that contains mutation His274Tyr and was identified as an oseltamivir-resistant strain[1],[17].

The reliability of homology model of N1 was identified by Ramachandran plot[11],[15]. It was limited by an orange area, that had coordinate of secondary protein structure as maximum tolerance limit area of steric strain (. This area was an allowed region, if there was a protein out of the blue line, putting the amino acid in the outlier (dissallowed region). In the tolerance limit area of residues plot, glysine was not allowed because it did not have a side chain, making φ (phi) and ψ (psi) angles not limited. The number of residue plot besides glysine showed the quality of the protein structure. The results of the Ramachandran plot showed that only Lys331 was in the outlier. About 95.8% of the residues were in the core region of the Ramachandran map, so that the homology model was used for the docking process. Superimposition was used to find out the degree of similarity between template 3CKZ chain A (PDB code) and the homology model. RMSD was calculated to find out the structural similarity between these proteins. The RMSD value of superimposition was 0.07 Å, indicating that the homology model N1 was very close to the structural similarity and was therefore used for the docking process.

Fig. 1

Ramachandran plot of modeling structure of N1. It shows only one outlier Lys331 in dissallowed region.

Screening using molecular docking

A total of 1,232 oseltamivir modified ligands and oseltamivir were screened using molecular docking based on the lowest binding energy (ΔG)[28]. The docking of the 1,232 ligands resulted in 100 best ligands; then, repeated docking resulted in 20 best ligands. Screening was based on Lipinski's rule of five, which showed that the molecular weight for drug likeness must be under 500 g/mol[43],[44].

Screening of 20 best ligands resulted in 10 ligands. Those ligands were screened, resulting in three of the best ligands. They were selected based on the quantities of hydrogen to the catalytic site of neuraminidase. Most of the three best ligands have the O-glycoside and the apioside group as the functional group of modified oseltamivir[17],[28]. These ligands were AD3BF2D, CA1G3B and F1G4B. The apioside group greatly enhanced inhibitory activity[28] because it is similar to O-glicoside compounds that have inhibitory activity against neuraminidase, and especially hydrolize influenza virus sialidase[27].

In , the properties of three of the best ligands AD3BF2D, CA1G3B, and F1G4B and oseltamivir as standard ligand showed that these ligands had the lowest binding energy than oseltamivir as standards. The level of the lowest ΔG value was AD3BF2D, CA1G3B and F1G4B ligands with ΔG  =  −7.8885, −7.6293, and −7.5637 kcal/mol, respectively. These values were better than oseltamivir as the standard ligand with ΔG  =  −4.2841 kcal/mol. pKi values of these ligands were 15.636, 12.689, and 13.630 μM, respectively. These values were better than pKi values of oseltamivir with 10.745 μM. These pKi values of AD3BF2D, CA1G3B, and F1G4B ligands indicated that these ligands were effective and had better affinity to strongly interact with neuraminidase than oseltamivir. The molecular weights of these ligands were under 500 g/mol. AD3BF2D, CA1G3B and F1G4B ligands had a number of bindings with the catalytic site (10, 10 and 5 bindings, respectively) (). These ligand bindings in the catalytic site were chosen for molecular dynamics simulation as oseltamvir had fewer bindings with the catalytic site[28],[45].

Table 1

The properties of the three best ligands and standard ligand

Ligand	AD3BF2D	CA1G3B	F1G4B	OTV (standard)
ΔG (kcal/mol)	−7.8885	−7.6293	−7.5637	−4.2841
pKi (μmol/L)	15.636	12.689	13.630	10.745
Molecular weight	489.542	462.496	481.475	313.418
H Donor	6	3	8	4
H Acceptor	7	11	3	1
Binding to the catalytic site	10	10	5	0

Hydrogen bonds between amino acid residues of neuraminidase with AD3BF2D, CA1G3B and F1G4B ligands and oseltamivir are shown in . AD3BF2D, CA1G3B and F1G4B ligands had more binding affinity with the catalytic site than oseltamivir. AD3BF2D ligand had hydrogen bonds with Arg118, two bindings with Glu278 and Arg 293, two bindings with Arg368 and four bindings with Tyr402. CA1G3B ligand had hydrogen bonds with Arg118 and Asp151, three bindings with Arg293, two bindings with Arg368 and three bindings with Tyr402. F1G4B ligand has hydrogen bonds with Asp151, two bindings with Glu278 and two bindings with Tyr402. However, oseltamivir did not have hydrogen bond for interacting with the catalytic site of neuraminidase. The quantities of the interacting hydrogen bonds with the catalytic site of neuraminidase indicated that these ligands have the ability to inhibit neuraminidase. Ligands AD3BF2D, CA1G3B and F1G4B had many hydrogen bond interactions with the catalytic site of neuraminidase. The docking poses in N1 of these three ligands are shown in . These hydrogen bonds were chosen as the basis of molecular dynamics simulation.

Table 2

Interactions of three of the best ligands and standard ligands with neuraminidase after docking

Ligand	Interaction to the neuraminidase
OTV (standard)	Glu119, Glu119, Arg156, Trp179, Glu228
AD3BF2D	Arg118, Glu119, Glu228, Glu278, Glu278, Arg293, Arg368, Arg368, Tyr402, Tyr402, Tyr402, Tyr402
CA1G3B	Arg118, Glu119, Asp151, Arg152, Arg152, Ser247, Arg293, Arg293, Arg293, Arg368, Arg368, Tyr402, Tyr402, Tyr402
F1G4B	Glu119, Asp151, Arg152, Ser247, Glu277, Glu277, Glu278, Glu278, Asn295, Tyr402, Tyr402

Red text: residues in catalytic site; black text: residues in framework.

The evolution time of the hydrogen bonds from the inhibitor-enzyme complex provides an approach to evaluate the convergence of the dynamical properties of the system[25]. The interaction of the best ligands, AD3BF2D, CA1G3B and F1G4B with the enzyme was evaluated in the hydrated state using molecular dynamics simulation at two different temperatures (300 K and 312 K). The simulation at temperature of 300 K was chosen to evaluate the interaction of the complex enzyme-ligand at the room temperature.

The simulation at temperature of 312 K was chosen to evaluate the interaction of the complex enzyme-ligand during the fever. As shown in , AD3BF2D has good interaction than the other ligands, CA1G3B and F1G4B. Compared with oseltamivir as the standard ligand ( and ), AD3BF2D had better interaction in two simulations at temperature of 300 K and 312 K. The interaction of AD3BF2D to form hydrogen bonds was much more stable in molecular docking, until the end of the simulations. However, the interaction between molecular docking and molecular dynamics simulation decreased from 10 hydrogen bonds to 1 and 3 hydrogen bonds at the end of simulations of temperatures of 300 K and 312 K, respectively.

Table 3

Interaction of ligands

Ligand	After docking	The end of simulations at temperature of 300K	The end of simulations at temperature of 312K
OTV (standard)	Glu119, Glu119, Arg156, Trp179, Glu228	Glu277, Glu278, Trp296,	Phe121, Arg225
AD3BF2D	Arg118, Glu119, Glu228, Glu278, Glu278, Arg293, Arg368, Arg368, Tyr402, Tyr402, Tyr402, Tyr402	Arg225, Glu228, Ser229, Ser229, Glu277, Glu278	Glu278, Arg293, Arg293,
CA1G3B	Arg118, Glu119, Asp151, Arg152, Arg152, Ser247, Arg293, Arg293, Arg293, Arg368, Arg368, Tyr402, Tyr402, Tyr402,	Glu119, Glu228, Glu277, Glu278, Arg293, Lys432,	Arg118, Glu228, Glu277, Glu278
F1G4B	Glu119, Asp151, , Arg152, Ser247, Glu277, Glu277, Glu278, Glu278, Asn295, Tyr402, Tyr402	Arg118, Glu119, Ala138, Glu277, Glu278, His297,	Ser180, His185, Glu230,

Red text: residues in catalytic site; black text: residues in framework.

Fig. 2

The docking processes of the three best ligands in N1. A: OTV, B: AD3BF2D, C: CA1G3B, and D: F1G4B.

Fig. 3

Interaction to form hydrogen bonds of OTV ligand at the end simulations of temperature of 300 K.

The stability of hydrogen bonds occured between Glu278 as one of the catalytic residues of neuraminidase and NH3 as the functional group of AD3BF2D ligand during simulations at temperature of 300 K (). The stability of hydrogen bonds also occured between Glu278 and NH3 functional group of AD3BF2D during simulations at 312 K. The stability of hydrogen bonds was also increased at simulations at 312 K between Arg293 and alcohol of the apioside group as the functional group of AD3BF2D ligand; and this interaction with Arg293 as catalytic residue of neuraminidase resulted in two hydrogen bonds (). If compared with the interaction that has 10 hydrogen bonds after docking, in the binding catalytic site, interaction after simulations could not form hydrogen bonds due to the mobility of amino acid residues in the catalytic site. It had different motions in solvent or hydrated conditions. These motions proved that there were dynamic enzymes due to the effect of solvent molecules. The presence of the different functional groups in each ligand had various bond effects that resulted in the different conformations. The stable interaction of AD3BF2D ligand to form hydrogen bonds during simulations at two different temperatures provided solid evidence that this ligand had affinity with neuraminidase. More complete ligand interactions data are in the supplementary material (http://www.bioinf.uni-leipzig.de/~arli/supplementary_material.pdf). As shown in , there were changes in conformation of the enzyme-AD3BF2D ligand complex, which occured during simulation, and showed the dynamic behaviour of the enzyme in the presence of solvent and ligand. AD3BF2D ligand remained in the binding pocket after docking, initialization stage and at the end of simulations at 300 K and 312 K. The binding pocket of this complex changed from the ‘closed’ form into the ‘open’ state at the end of simulation. Based on the simulation results, the changes in conformation of enzyme-AD3BF2D ligand complex during simulations caused apioside and ammonium functional groups to form more stable interactions with Glu278 and Arg293 compared with others. Those functional groups, namely pentyloxy, aldehyde, and carboxylate did not form hydrogen bonds with binding pocket. More conformation figures of other ligands are listed in the Supplementary Material). shows that the ligand-receptor interaction indeed reached its stability after 2,000–3,000 ps. This would suggest that the ligands, as drug candidates, are within the stable condition in the human body.

Fig. 4

Interaction to form hydrogen bonds of OTV ligand at the end simulations of temperature of 312 K.

Fig. 5

Interaction to form hydrogen bonds of AD3BF2D ligand at the end simulations of temperature of 300 K.

Fig. 6

Interaction to form hydrogen bonds of AD3BF2D ligand at the end simulations of temperature of 312 K.

Fig. 7

The RMSD Graphic of (A) CA1GB and Neuraminidase (NA) dynamics at 300 K.

The X axis represents the molecular dynamics time duration in picoseconds (ps), while the Y axis represents the RMSD value. The graphics shows that after 2,000 ps, the structure become stable. B: The RMSD Graphic of CA1GB and Neuraminidase (NA) dynamics at 312 K. The graphics shows that after 2,000 ps, the structure become stable. C: The RMSD Graphic of F1G4B and Neuraminidase (NA) dynamics at 300 K. The graphics shows that after 2,000 ps, the structure become stable. D: The RMSD Graphic of F1G4B and Neuraminidase (NA) dynamics at 312 K. The graphics shows that after 2,000 ps, the structure become stable. E: The RMSD Graphic of AD3BF and Neuraminidase (NA) dynamics at 300 K. The graphics shows that after 3,000 ps, the structure become stable. F: The RMSD Graphic of AD3BF2D and Neuraminidase (NA) dynamics at 312 K. The graphics shows that after 3,000 ps, the structure become stable.

Discussion

In this research, the contribution from functional groups AD3BF2D ligand (N-[(1S,6R)-5-amino-5-{[(2R,3S,4S)-3,4-dihydroxy-4-(hydroxymethyl) tetrahydrofuran-2-yl]oxy}-4- formylcyclohex-3-en-1-yl]acetamide-3-(1-ethylpropoxy)-1-cyclohexene-1-carboxylate), apioside and ammonium group provide better interaction to form hydrogen bonds as affinity to neuraminidase than oseltamivir due to the polarity of functional groups, so that these functional groups are able to form interaction with neuraminidase that have polar and hydrophilic residues in the catalytic site. The apioside group forms two hydrogen bonds with Arg293, and the ammonium group forms one hydrogen bond with Glu278 as a residue in the catalytic site of neuraminidase. This result of interaction of the apioside and ammonium group with the residue in the catalytic site of neuraminidase indicated that AD3BF2D ligand can inhibit neuraminidase and this ligand can be proposed as a candidate for neuraminidase inhibitor of influenza A virus subtype H1N1.

As an outlook of this study, we are aware that some research groups already produced RNA based drugs for interfering with influenza virus[46]-[48]. Moreover, a feasible and powerful transcriptomics based tool has been generated for viral bioinformatics research[49]-[52]. To this end, it would be interesting to investigate the transcriptomics properties of the H1N1 virus, as supplementary approach to the existing proteomics method.

In conclusion, we have built the latest N1 structure model by homology modeling, which has excellent reliability by using Ramachandran plot and template superimposition. A total of 1,232 oseltamivir modified ligands molecules have been designed and screened by molecular docking. Three of the best ligands, AD3BF2D, CA1G3B, and F1G4B were obtained with the lowest binding energy that contradicted to OTV as standard. These ligands' interactions were evaluated in hydrated state using molecular dynamics simulation. This work has provided AD3BF2D ligands that have stable interaction with Glu278 and Glu 278, Arg293 and Arg 293 during simulations at temperature of 300 K and 312 K, respectively. Based on these results, AD3BF2D ligand can be used as a candidate for neuraminidase inhibitor against influenza virus A subtype H1N1.

Acknowledgements

The authors would like to thank Hibah BOPTN PUPT Ditjen Dikti No: 0970/H2.R12/HKP.05.00/2014 and Directorate of Research and Community Engagement, University of Indonesia, for supporting this research. Usman Sumo Friend Tambunan supervised this research, Rizky Archintya Rachmania was working on the technical details and preparing the research report and Arli Aditya Parikesit was responsible for writing and revising the manuscript.