IN SILICO DESIGN OF THE M2 PROTON CHANNEL INHIBITORS OF H1N1 VIRUS

M2 proton channel of H1N1 Influenza A virus is the target protein anti-flu drugs amantadine and rimantadine. However, the two once powerful adamant ane-based drugs lost their 90% bioactivity because of mutations of virus in recent twenty years. The resi stance of influenza A virus to amantadine need to d evelop more effective adamantane-based drugs. Several rese arch by molecular docking method have been conducted to design and discover ligands which beco m potential inhibitors for the M2 channel protein of influenza virus in order to inhibit the replication of influenza virus. In this research are studied a n evaluated the interaction of ligands towards the pr otein in the hydrated state using molecular dynamic s simulations at two different temperatures. Analysis of ligand interaction yields that AM-L6-R6 ligand has best affinity towards the protein than the T-R6-L6, T-L6-R12 and the standard ligand. It is shown by t he ligand interaction on the enzyme active site which remains to be formed during the simulation performe d. At the end of simulation temperature of 300 K, AM-L 6 R6 ligand has a residue contact with the Arg45 an d formed hydrogen bond with Asp44. Then at the end of simulation temperature of 312 K, AM-L6-R6 ligands also could form a hydrogen bond with Asp44. Conform ational changes of protein which occur during simulation showed the dynamicization of an protein in the presence of solvent and inhibitor.


INTRODUCTION
Recently, the outbreak of H1N1 influenza A virus is a pandemic of a new strain of influenza virus identified in April 2009 (Trifonov et al., 2009), commonly to as "swine flu". Within only four months, the pandemic has caused many deaths from the first detected country Mexico to almost all countries of the world (http://www.who.int/csr/disease/swineflu/). The H1N1 influenza virus is quite familiar to us because it had caused the 1918-1919 Spain pandemic that had infected 5% of the world population and resulted in 20-50 million deaths worldwide (Trifonov et al., 2009). In July 2009 the World Health Organization (WHO) enhanced the warning phase 6, meaning that the spread of H1N1 influenza virus has become a serious global pandemic. The even worse news is that cases were reported that several strains of H1N1 influenza a viruses were resistant to oseltamivir (Tamiflu) (Du et al., 2009).
The M2 channel of the influenza A virus is a homotetrameric protein that contains three different kinds of Transmembrane (TM) four-helix channel with 97 residues per subunit, each of which comprises an intracellular C-terminal domain with a length of 54 residues, a transmembrane domain with a length of 19 residues and an extracellular N-terminal domain with a length of 24 residues (Holsinger and Lamb, 1991;Sugrue and Hay, 1991;Huang et al., 2008;Du et al., 2009;Intharathep et al., 2008). This protein has a major function as a proton-selective channel that is controlled Science Publications OJBS by the endosomal pH values (Shimbo et al., 1996;Huang et al., 2008;Du et al., 2009;Intharathep et al., 2008). In conducting endosomal protons into the virion, the M2 channel is considered to play a crucial role in the viral life cycle (Bright et al., 2006;Wang et al., 2009). Once protons flux into the virion, the uncoating process of the viral nucleic acid in endosomes will be activated, which is essential for the virion to invade the host cells. The pH regulation of the influenza virus M2 channel is believed to have strong relationship with the TM fourhelix bundle (Chou et al., 1988; that is experimentally demonstrated to have proton conductive activity (Duff and Ashley, 1992;Chen et al., 2007). Additionally, this channel has another function to preserve the pre-fusion state of HA during viral maturation Chen et al., 2007;Grambas and Hay, 1992;Takeuchi and Lamb, 1994).
Influenza A virus has the ability to undergo changes by the mechanisms of antigenic drift and shift, resulting new evolving virus strains, which may be extremely toxic and drug resistant (Du et al., 2010;Deyde et al., 2007;Smith et al., 2002). Given that influenza shift may occur every 20-30 years, the danger of future influenza a pandemics highlights the need to develop more effective drugs. The threat of an impending influenza pandemic, possibly through the mutations of the present avian strain H5N1 or swine strain H1N1, has triggered a global effort to develop more effective antivirus drugs. However, during the past several decades many efforts in developing anti influenza drugs have almost been futile due to the rapid mutations of the influenza virus, resulting in the persistent resistance to the exixting drugs (Du et al., 2010).
Although vaccination is the ideal way to prevent influenza virus infection, the preparation of a new vaccine requires more than 6 months (Couzin-Frankel, 2009). Thus, antiviral drugs are the most effective for short-term defense against influenza. The only known anti-influenza A drugs are M2 (matrix-2-protein) inhibitors (amantadine and its derivative rimantadine) and NA (neuraminidase) inhibitors (zanamivir and oseltamivir) (Hayden and Hay, 1992;Schmidtke et al., 2006).The mechanism of M2 inhibitors is to block the ion channel activity of the M2 protein of most influenza A viruses. This action inhibits viral replication by blocking proton flow. The amino group in amantadine is likely the pharmacophore and is necessary to block proton transport (Hatta et al., 2004;Wang et al., 1993;Betakova, 2007;Cady and Hong, 2008;Tran et al., 2011). Both amantadine and rimantadine are limited in their use in the U.S because of the rapid development of resistance. In addition, there is growing concern that antineuraminidase-resistant viruses may emerge if these drugs are widely used (Lipatov et al., 2004). Thus, there is an urgent need to discover new types of M2 inhibitors for the development of new antiinfluenza drugs. The present study was initiated in an attempt to solve drug resistant problem, to design and screen a small primary amine library of scaffold-hops based on amantadine to generate new hits in the M2 inhibitors class by conducting molecular docking and dynamics method using MOE software (Vilar et al., 2008).

MATERIALS AND METHODS
The following steps were conducted using Microsoft Windows 7 based PC. Search and choose the sequences (Veits et al., 2012). M2 proton channel sequences were downloaded from Influenza Virus Resource database of the National Center for Biotechnology Information (NCBI) (http://www.ncbi.nlm.nih.gov.html).

Multiple Sequence Alignment and Homology Modeling
This step was conducted by using ClustalW online program (www.ebi.ac.uk/Tools/clustralw2/index.html). Moreover, Homology modeling was performed using the SWISS MODEL which can be accessed through http://swissmodel.expasy.org/SWISS-MODEL.html. The SWISS MODEL service was utilized to build 3D models for H1N1 virus protein, by finding the exact templates (Abdussamad and Aris-Brosou, 2011). The PDB files were exposed by using Molecular Operating Environment (MOE) (Vilar et al., 2008).

Molecular Docking
The found 3D structure was docked with its ligand by using molecular docking software. Before the docking, preparation steps must be done. This was done by removing water molecules, addition of hydrogen atoms and charges. Further minimization was done using MMFF94x (Vilar et al., 2008). The utilized parameters for analyzing the complex between protein and ligand are ∆G° binding and inhibition constant.

Molecular Dynamics
The simulation on protein/ligand complex was done after molecular docking steps with MOE 2008.10. Before molecular dynamics was computed preparation steps were done for molecular docking, followed by inserting the ligand in order to form the protein-ligand complex. Then, the complex was minimized by force field MMFF94x and solvable in Born form. The parameters were utilized in accordance with MOE default, which is ensemble NVT (N,

Science Publications
OJBS total atom; V, Volume; T, Temperature) by using the NPA algorithm. We employed two different temperatures level in this simulation, 300 K as default MOE temperature and 312 K for fever temperature. The system was heated up to 312 K through a constant-volume MD simulation lasting 20 ps. The enzyme-inhibitor binding interactions and stability of enzyme-ligand complex were calculated after 5 ns of simulation for all ligands (Tambunan et al., 2010).

Multiple Sequence Alignment
Sequence alignments performed using ClustalW2 program that can be accessed through European Bioinformatics Institute (EBI) (http://www.ebi.ac.uk/Tools/clustalw2/index.html). The downloaded M2 channel FASTA sequences were loaded into ClustalW2 program for multiple alignment process. Sequence alignments were used to see differences of protein sequences and to determine a single representative of M2 channel that will be targeted.

Homology Modeling
The sequence alignment was performed to determine template for docking simulation. Swiss-model (http://swissmodel.expasy.org/) was employed for template identification for M2 channel (Du et al., 2010;Arnold et al., 2006). The process was conducted by giving sequence file input for browsing the crystal structure in FASTA format. The chosen mode is automated. Based on template identification result, crystal structure of the M2 channel (PDB entry code: 2rlfD) were utilized.

Molecular Docking
The protein structure was prepared for docking (Pyrkov et al., 2010). This included the addition of missing hydrogens (Gordon et al., 2005). During docking, series of poses (ligand-protein complexes) were generated for each molecule. The quality of each pose was further assessed by the London dG (LdG) Scoring Function (SF) which estimated the binding free energy of the ligand and a set of the highest score poses were chosen for each molecule (pose) docked (Mazur et al., 2009;Vilar et al., 2008). In this study, MMFF94 and MMFF94x were used for forcefield minimization. MMFF94x was reported as the efficient forcefield for minimizing ligand-protein complexes (Mazur et al., 2009;Kaminski and Jorgensen, 1996;Halgren, 1990). Triangle Matcher (Wildman and Crippen, 2001;Cook et al., 2010) was applied to orient the ligands in the active site based on charge groups and spatial fit. Triangle Matcher performs a random walk with the ligand in the active site to define the optimal binding orientations. Based on docking simulation result, there are three ligand, namely AML6R6, TR6L6 and TL6R12, could be proposed as a potential inhibitor to the interaction of M2 channel in A/2009(H1N1) virus (Fig. 2).

OJBS
In Table 2, it shows the calculated free binding energy (∆Gbind) and hydrogen bond of flexible-ligand docking simulation from the best ligand (Frimurer et al., 2003). Free energy of binding for all the designed ligands were lower than the comparative ligand (amantadine and rimantadine) (Du et al., 2010). The free energy of AML6R6, TR6L6, TL6R12, rimantadine and amantadine ligands were negative value, respectivelly, -6. 6766,-5.7709, -4.9738, -4.1198 and -3.9574 kcal/mol. The negative and low value of ∆G bind indicated the strong favorable bond between enzyme and ligand. In Table 2, we also presented hydrogen bond of AML6R6, TR6L6, TL6R12, rimantadine and amantadine ligands. Bold letters indicate residues of the active site of the M2 channel. The results showed that AML6R6, TR6L6 and TL6R12 binds at channel lock Asp44. Illustrated in Fig. 3-5, ligand T-R6-L6 and L6-R12-T form two hydrogen bonds with functioanal residues Asp44. The Asp44 residue is not only the channel lock, but also the protein exit of the M2 channel. Therefore, drug binding at this position not only can hold the channel in the closed conformation but also can simply block the proton exits (Pielak et al., 2009). While the standard ligand amantadine and rimantadine do not form hydrogen bonds with the active site of the M2 channel. Interactions that occur will be validated with the state of the hydrated molecules using molecular dynamics simulation method (Du et al., 2010;Kumar, 2011;Jing et al., 2008;Zhou and Schulten, 1996).

Molecular Dynamics Simulation
Molecular dynamics simulations are one of the most versatile and widely applied computational techniques for the study of biological macromolecules. They are very valuable for understanding the dynamic behavior of proteins at different timescales, from fast internal motions to slow conformational changes or even protein folding processes (Snow et al., 2005;Alonso et al., 2006). It is also possible to study the effect of explicit solvent molecules on protein structure.
In Table 3, we presented comparison hydrogen bond of AML6R6, TR6L6, TL6R12, rimantadine and amantadine ligands after molecular docking and dynamics simulation. Molecular docking process has resulted all shown ligand to form hydrogen bonds with the functional residues Asp44 (Du et al., 2010). After the simulation conducted, the three ligands are still forming hydrogen bonds with Asp44. AML6R6 and TL6R12 ligands form two hydrogen bonds with Asp44 residue after molecular dynamics simulations (Huang et al., 2008). However TR6L6 ligand losing one hydrogen bond with Asp44 after molecular dynamics simulation process.  Residue contact is interaction of M2 channel residue between ligand (inhibitor) (Du et al., 2010). Interaction between complex ligand-M2 channels is formed with hydrogen bond, but also it occurs non-covalent interactions which influence the stability of proteins (Ball and Maechling, 2009). Hydrophobic interactions that cause the protein to maintain its integration because of folding and aromatic interactions that affect the stability of tertiary structure. In Table 3 shown AML6R6, TR6L6 and TL6R12 have more residue contacts (italic letter) with M2 channel after molecular dynamics simulations. TL6R12 ligand has more residue contact (Arg45, Phe47, Phe47, Phe48) than two other ligands.

DISCUSSION
The development of molecular docking and dynamics method has eventually making the drug design workable. The advances of the utilized software are sufficient for generating the necessary ligands and then redirect them for the molecular docking and dynamics computation. Moreover, the pipeline is sufficient for handling this huge amount of data. The time and space complexity of this pipeline are handled by the efficacy of the available graphic subsystem of existing hardware.
Our pipeline could cope with more than 1400 lead compounds, which are much more numerous than what we have before. The screening pipeline was resulted in only three potential inhibitors that perform much better than the standards. Here, the scaling up of our developed pipeline to cope with bigger amount of data seemed to be feasible, provided that the existing hardware and software architecture will support it.
Hydrogen bond, as one of the important chemical bond, is definitely a focal point of this study due its important role for maintaining the interaction between ligand and the protein. The calculation of binding free energy and the related binding sites have eventually exposed important clues on the reactivity of the drugs candidate.
The determination of the molecular dynamics at 300 and 312 K are indeed useful for simulating the normal and fever human body temperature. The generated data has already sufficient for predicting the stability of the drugs candidate in both temperatures. In this end, a better predictability of the drug mechanism could be achieved by the improved simulation of the temperature.

CONCLUSION
There are 1447 amantadine-based drugs already designed in this study. The docking of these inhibitors to determine their interaction activities with M2 channel was carried out. The free energy binding (∆G bind ) values for all the inhibitors were small, indicating the molecules bound considerably well to the binding site. From docking process there are three ligands that have smaller free binding energy (∆G bind ) than ligand standard (amantadine and rimantadine) and the best interaction with M2 channel, there are AML6R6, TR6L6 and T-L6-L6-R12. Docking result shown that all ligands form hydrogen bond with functional residue Asp44. After dynamics simulation at 300 K and 312 K, all ligands still form hydrogen bond with Asp44 residue. To ensure the absorption, distribution, metabolism, excretion and toxicity, ADMET calculations on the ligands are recommended. Thus, from this study, the AML6R6, TR6L6 and TL6R12 could be proposed as a potential inhibitor to inhibit interaction of M2 channel in H1N1 virus.

Competing Interest
The researchers declare that they have no competing interests.

ACKNOWLEDGEMENT
The research was supported by Grant Hibah Bersaing Pasca Sarjana, Indonesian Ministry of Education (Contract Number: 2848/H2.R12/PPM.00 Penelitian/2012). The researchers are grateful to Asif M. Khan, Perdana University, Selangor Darul Ehsan, Malaysia for proof reading the manuscript. USFT supervised this research, MR was working on the technical details and AAP was preparing the English manuscript and re-verified the data.