key: cord-0001772-3gz4dyie authors: Tambunan, Usman Sumo Friend; Rachmania, Rizky Archintya; Parikesit, Arli Aditya title: In silico modification of oseltamivir as neuraminidase inhibitor of influenza A virus subtype H1N1 date: 2014-12-12 journal: J Biomed Res DOI: 10.7555/jbr.29.20130024 sha: 48f7616ac0756fa4bcdb79fef569643b929b108e doc_id: 1772 cord_uid: 3gz4dyie 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. 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] [2] [3] [4] [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] [10] [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 proteinligand interactions with atomic resolutions based on classical mechanics [22] [23] [24] [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. 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 rootmean-square deviation (RMSD) was calculated to identify structural similarity between template model mutated with 3D structure from homology modeling. 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 (DG binding ) from pose of ligand to enzyme was based on the calculation: Where, c 5 rotation entropy rate and translation that can be obtained or released, E flex 5 flexibility of ligand energy, f HB 5 imperfection size geometry of hydrogen bonds, c HB 5 energy of an ideal hydrogen bond, f M 5 imperfection size geometry of metal ligations, c M 5 ideal metal ligation energy and D i 5 atom desolvation energy [31, [34] [35] [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 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 6A. Dynamics simulation stage involved initialization for 40 ps, main simulations as equilibrium and production stage for produced trajectory. Molecular dynamics simulation of complex enzymeligand 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 (E sol ) to be calculated in potential energy molecular system function from atom coordinate: Where, E str 5 stretching energy, E ang 5 angular energy, E stb 5 stretching and bending energy, E oop 5 operating energy, E tor 5 energy torsional, E vdw 5 intermolecules van der Waals energy, E ele 5 electrostatic energy, E sol 5 solvation energy, dan E res 5 residual moving energy so that we can see clearly Esol in the equation [39] [40] [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. 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 (Fig. 1) . 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 Q (phi) and y (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 A, indicating that the homology model N1 was very close to the structural similarity and was therefore used for the docking process. A total of 1,232 oseltamivir modified ligands and oseltamivir were screened using molecular docking based on the lowest binding energy (DG) [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 9 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 Table 1 , 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 DG value was AD3BF2D, CA1G3B and F1G4B ligands with DG 5 -7.8885, -7.6293, and -7.5637 kcal/mol, respectively. These values were better than oseltamivir as the standard ligand with DG 5 -4.2841 kcal/mol. pKi values of these ligands were 15.636, 12.689, and 13.630 mM, respectively. These values were better than pKi values of oseltamivir with 10.745 mM. 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) ( Table 1) . These ligand bindings in the catalytic site were chosen for molecular dynamics simulation as oseltamvir had fewer bindings with the catalytic site [28, 45] . Hydrogen bonds between amino acid residues of neuraminidase with AD3BF2D, CA1G3B and F1G4B ligands and oseltamivir are shown in Table 2 . 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 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 Fig. 1 . These hydrogen bonds were chosen as the basis of molecular dynamics simulation. 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 Table 3 , AD3BF2D has good interaction than the other ligands, CA1G3B and F1G4B. Compared with oseltamivir as the standard ligand ( Fig. 2 and 3) , 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. The stability of hydrogen bonds occured between Glu278 as one of the catalytic residues of neuraminidase and NH 3 as the functional group of AD3BF2D ligand during simulations at temperature of 300 K (Fig. 4) . 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 (Fig. 5) . 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 Fig. 6 , 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). Fig. 7A -F 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. 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. 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] [47] [48] . Moreover, a feasible and powerful transcriptomics based tool has been generated for viral bioinformatics research [49] [50] [51] [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 9 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. Insights from investigating the interaction of oseltamivir (Tamiflu) with neuraminidase of the 2009 H1N1 swine flu virus Three-dimensional structure of the complex of 4-guanidino-Neu5Ac2en and influenza virus neuraminidase Mechanism by which mutations at his274 alter sensitivity of influenza a virus n1 neuraminidase to oseltamivir carboxylate and zanamivir Drug screening for influenza neuraminidase inhibitors Analogs of zanamivir with modified C4-substituents as the inhibitors against the group-1 neuraminidases of influenza viruses Susceptibility of antiviral drugs against 2009 influenza A (H1N1) virus Mutations of neuraminidase implicated in neuraminidase inhibitors resistance Global transmission of oseltamivir-resistant influenza Design of Linear Peptide as Neuraminidase Inhibitor Influenza A Virus Base on Molecular Docking Simulation Drug design for Influenza A virus subtype H1N1 Insights from modeling the 3D structure of H5N1 influenza virus neuraminidase and its binding interactions with ligands Emergence and spread of oseltamivir-resistant A(H1N1) influenza viruses in Oceania, South East Asia and South Africa Widespread oseltamivir resistance in influenza A viruses (H1N1), South Africa Novel inhibitor design for hemagglutinin against H1N1 influenza virus by core hopping method Analogue inhibitors by modifying oseltamivir based on the crystal neuraminidase structure for treating drug-resistant H5N1 virus In Vitro Antiviral Activity of AG7088, a Potent Inhibitor of Human Rhinovirus 3C Protease In silico design of cyclic peptides as influenza virus, a subtype H1N1 neuraminidase inhibitor Oseltamivir (Tamiflu) and its potential for use in the event of an influenza pandemic Importance of neuraminidase active-site residues to the neuraminidase inhibitor resistance of influenza viruses Structural bioinformatics and its impact to biomedical science Binding mechanism of coronavirus main proteinase with ligands and its implication to drug design against SARS Molecular dynamics: survey of methods for simulating the activity of proteins Biomolecular modeling: Goals, problems, perspectives Structure and localisation of drug binding sites on neurotransmitter transporters Three new powerful oseltamivir derivatives for inhibiting the neuraminidase of influenza virus Prediction of antimicrobial peptides based on sequence alignment and feature selection methods An O-glycoside of sialic acid derivative that inhibits both hemagglutinin and sialidase activities of influenza viruses Inhibition of neuraminidase activity by polyphenol compounds isolated from the roots of Glycyrrhiza uralensis ACD Labs/LogP dB 3.5 and ChemSketch 3.5 Screening of Bioactive Compounds from Olea Europaea as Glutamine Synthase A Inhibitor of Bacterial Meningitis haemophilus Infulenza Type B Through Molecular Docking Simulation In silico Modification of (1R, 2R, 3R, 5S)-(-)-Isopinocampheylamine as Inhibitors of M2 Proton Channel in Influenza A Virus Subtype H1N1, using the Molecular Docking Approach Medicinal chemistry and the molecular operating environment (MOE): application of QSAR and molecular docking to drug discovery Bioactive Compounds Screening from Zingiberaceae Family as Influenza A/Swine Flu Virus Neuraminidase Inhibitor through Docking Approach Solvated interaction energy (SIE) for scoring protein-ligand binding affinities. 1. Exploring the parameter space Generalized Born Model: Analysis, Refinement, and Applications to Proteins A computational approach to evaluate the androgenic affinity of iprodione, procymidone, vinclozolin and their metabolites Molecular dynamics approaches estimate the binding energy of HIV-1 integrase inhibitors and correlate with in vitro activity Computational design of disulfide cyclic peptide as potential inhibitor of complex NS2B-NS3 dengue virus protease Molecular Dynamics Simulation of DENV RNA-Dependent RNA-Polymerase with Potential Inhibitor of Disulfide Cyclic Peptide Simulations of solvation dynamics in simple polar solvents Silico Detection, Drug Design and Prevention Agent Development. In: Rajkumar R, ed. Topics on Cervical Cancer with an Advocacy for Prevention The high-affinity binding site for tricyclic antidepressants resides in the outer vestibule of the serotonin transporter Silico Analysis of Envelope Dengue Virus-2 and Envelope Dengue Virus-3 Protein as the Backbone of Dengue Virus Tetravalent Vaccine by Using Homology Modeling Method Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings In silico modification of suberoylanilide hydroxamic acid (SAHA) as potential inhibitor for class II histone deacetylase (HDAC) RNA-based therapeutics: ready for delivery? The promises and pitfalls of RNAinterference-based therapeutics Therapeutic Potential of RNA Interference Automatic detection of conserved RNA structure elements in complete RNA virus genomes Conserved RNA secondary structures in Flaviviridae genomes Automatic detection of conserved base pairing patterns in RNA virus genomes Secondary structure prediction for aligned RNA sequences 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.