Discovery of Potential M2 Channel Inhibitors Based on the Amantadine Scaffold via Virtual Screening and Pharmacophore Modeling

The M2 channel protein on the influenza A virus membrane has become the main target of the anti-flu drugs amantadine and rimantadine. The structure of the M2 channel proteins of the H3N2 (PDB code 2RLF) and 2009-H1N1 (Genbank accession number GQ385383) viruses may help researchers to solve the drug-resistant problem of these two adamantane-based drugs and develop more powerful new drugs against influenza A virus. In the present study, we searched for new M2 channel inhibitors through a combination of different computational methodologies, including virtual screening with docking and pharmacophore modeling. Virtual screening was performed to calculate the free energies of binding between receptor M2 channel proteins and 200 new designed ligands. After that, pharmacophore analysis was used to identify the important M2 protein-inhibitor interactions and common features of top binding compounds with M2 channel proteins. Finally, the two most potential compounds were determined as novel leads to inhibit M2 channel proteins in both H3N2 and 2009-H1N1 influenza A virus.


Introduction
A functional M2 channel protein is essential for the release of the flu virus' genetic material inside an infected cell [1]. The M2 channel protein, which is a pH-sensitive proton channel, also plays a key role in virus replication and regulates virus morphology [2,3]. The two adamantane derivative-based drugs amantadine and rimantadine (Figure 1), which have been used as the first-choice antiviral drugs against community outbreaks, were the first antivirals licensed for effective against influenza A viruses [4]. However, since 2003, the frequency of amantadine resistance has increased significantly, from less than 5% to greater than 90% of isolated influenza A virus [5]. There is therefore a great urgent need to discover new types of M2 inhibitors for the development of new anti-influenza drugs due to the new mutations on M2 channel protein and the drug-resistance of amantadine and rimantadine [6]. The new drugs directed against M2 channel proteins should be more universal and effective than the previous ones.  Why have amantadine and rimantadine become resistant to influenza A virus recently? Site-directed mutagenesis and molecular dynamics simulations have been carried out to investigate the amantadine resistance in the trans-membrane domain of the M2 channel protein [7,8]. According to statistical data on resistant mutants, 70% to 80% of substitutions occur at position Ser31, and around 10% occur at positions Val27 and Ala30 in vitro and in clinical samples [7]. To solve the drug-resistance problem, a reliable molecular structure of M2 channel protein is a high priority for designing new drugs [10]. The M2 channel protein structures obtained experimentally in previous studies [11][12][13] have thus become the main targets for scientists and pharmacologists to find drugs against influenza A virus using structure-based drug design approaches [9,10]. One of them is the high-resolution nuclear magnetic resonance (NMR) spectroscopy structure by Schnell and Chou with the Protein Data Bank (PDB) code entry of 2RLF [14] that has successfully provided a full-length structure of H3N2 M2 channel protein. The 3D 2009-H1N1 M2 channel protein [15] built from sequence with the Genbank accession number of GQ385303 was also used in this current research.
In previous studies, drug binding affinities which modified the functional groups on amantadine did not reveal any details on how the ligands actually bind at the molecular level [16][17][18]. This research aims to search for drug candidates that are effective against the resistant strains of influenza A viruses and shed light on several important insight top hit M2 protein-inhibitor interactions. In this study, 200 drug candidates were designed by modifying or adding more functional groups to the amantadine scaffold and then used for virtual screening process [19]. After that, top 10 binding compounds were selected for further studied in pharmacophore analysis.

Binding Site Identification
Two possible binding sites for the M2 channel protein of influenza found in experimental studies are the drug binding locations [20]. The molecular docking results on both amantadine and rimantadine positioned inside and outside the M2 channel proteins partially supported the actual binding site location. The free energy of binding of amantadine and rimantadine inside the channel is generally lower than the binding outside the M2 channel proteins (cf. Table 1). Comparing the two cases, the affinity of the inside channel binding is higher than that in the outside of the channel. The ligands binding inside the H3N2 M2 channel protein bind through a hydrogen bond with Ala30 (cf. Figure 2) but no formation of hydrogen bond with 2009-H1N1 M2 channel protein was observed for either amantadine or rimantadine. Meanwhile, the outside channel binding of ligands in both M2 channel proteins was largely unchanged in terms of the free binding energy and the location of the binding sites. This indicated that the binding between the ligands and the M2 protein structure is energetically more favourable and stable with the ligands (either amantadine or rimantadine) inside the channel than outside the M2 channel proteins.
Since the exact location of the functional adamantine binding site has been a source of controversy [20], for better design of novel M2 inhibitors and further comparison analysis, two binding positions were taken into consideration in this study. (A) Two hydrogen bonds were found between the carboxyl group of Ala30 and the amino group of amantadine inside the H3N2 M2 channel protein, with bond lengths of 2.16 and 2.65 Å, respectively. This inside binding mechanism keeps amantadine in a stable channel-bound conformation, with consistent stoichiometry and agrees well with resistance mutants; (B) two hydrogen bonds (2.06-2.09 Å) were formed between the carboxyl group of Ala30 and the amino group of rimantadine. Hydrophobic interactions between methyl group and Ala30, Ile33 and Gly34 inside H3N2 M2 channel protein were observed.

Docking Results and Ranking Top Hit Compounds Binding Inside and Outside M2 Channel Proteins
The 200 compounds were divided into 7 groups based on the structure similarity. Their free energy of binding of 200 compounds to both H3N2 and 2009-H1N1 M2 channel proteins at two positions, namely inside and outside the M2 channel, is shown in the Table 2.                The selection of top binding compounds was mostly based on choosing the common compounds in both positions which have the lowest free energy of binding with the 2009-H1N1 M2 channel protein.
In comparison with amantadine and rimantadine, the top 10 compounds were ranked based on the lowest free energy of binding obtained from docking which indicates a better binding affinity. Two groups were selected from the lowest energy of binding from virtual screening inside and outside the M2 channel proteins, respectively. Table 3 shows that the affinity of M2 drug candidates binding inside the channel is much higher than binding outside the M2 channel proteins.  The group inside the channel generally binds with high affinity through hydrophobic interactions inside the M2 channel proteins, specifically to hydrophobic residues Ile33, Val27, Ala30. In particular, for compound I9, a high possibility of forming hydrogen bonds was predicted for the amino group of I9 and two residues Gly34 or Ile34, whereas the methyl groups interact with the hydrophobic part of the M2 channel proteins (Ala30 and Ile33). Subsequently, I9, with high affinity, was selected for further pharmacophore analysis of the top hits' M2 protein-inhibitor interactions. Besides, the group binding on the outside channel basically binds with lower affinity through hydrogen bonds with Asp44 in every compound in either of the M2 protein types, and some hydrophobic interactions with Arg45, Leu46, Phe47 and Phe48. Interestingly, both compounds I5 and I9 are identical to O9 and O3 respectively (cf. Table 2). Therefore, these two compounds were selected to be top hits for either the inside or outside binding M2 channel proteins for further pharmacophore analysis of specific drug target-protein interactions. Moreover, the free energy of binding of the top 10 compounds binding inside and outside the M2 channel proteins was compared in detail (cf. Figure 3).

Pharmacophore Analysis for Top hits M2 Channel-Inhibitor Interactions
LigandScout presents the interactions between M2 channel proteins and ligands, as well as with some excluded volume spheres corresponding to their 3D complex structures. In terms of potential drug-target binding in this study, the top two selected hit compounds were visualized clearly for the possible interactions with critical residues in the M2 channel proteins. The generated four pharmacophore models of two highly potent binding compounds both inside and outside the M2 channel proteins with their geometrical constraints and active sites were represented in Table 4. Table 4. Structure-based hypotheses were superimposed on the active site of 3D structure of M2 channel proteins. Hydrogen Bond Acceptor (HBA) is shown as green vectored spheres, Hydrophobic (H) as light blue spheres (for interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) As shown in Table 4, for the I5 inside the M2 channel protein complex, the generated pharmacophore contains one hydrophobic (H) chemical feature which points towards Ala30 and Ile33. A four features hypothesis was generated from the I9 inside the M2 channel complex, which comprises one HBA towards Ile33 and Gly34, and three H groups mostly oriented towards Ala30 with 14 excluded volume spheres. A five features hypothesis was generated from the I5 outside the M2 channel complex, including one HBA toward Asp44 and four H groups pointed towards Leu40 and Leu43 with 11 excluded volume spheres. The I9 outside the M2 channel complex produced a four features hypothesis consisting of one HBA and three H groups with nine excluded volume spheres. The HBA group pointed towards the Asp44 while the H groups pointed towards Arg45 and Leu46, respectively. Comparing the above four pharmacophore models, the HBA and H groups from all the models were pointed towards Ala30, Ile33, Gly34 inside the M2 channel protein and Asp44, Leu40, Leu43, Arg 45 and Leu46 outside the M2 channel, which plays a major role in inhibiting M2 channel protein activity, respectively. The results suggested that these residues are extremely important binding sites for novel M2 channel inhibitors.

Common Pharmacophore Features of the Top Binding Compounds
The best hypothesis from Hip-Hop was chosen from seventeen hypotheses of the inside group and eighteen hypotheses from the outside group using Common Feature Pharmacophore Generation/Discovery Studio. As a result, the two representative pharmacophore features of the top 10 inside binding and outside binding compounds were identified (cf. Figure 4). For the group of compounds binding inside the M2 channel proteins, the generated pharmacophore contains one HBA and three H chemical features. Similarly, a three features hypothesis was generated from the group of outside binding compounds, which was composed of one HBA and two H groups. Hence, HBA and H features are considered as important chemical features to discover novel M2 channel inhibitors. Consistently with what predicted was by LigandScout, common features from Discovery Studio also indicate that all top binding compounds either inside or outside fit the common features having the same pharmacophore functional groups that can interact with critical residues in both the H3N2 and 2009-H1N1 M2 channel proteins.

Proteins and Inhibitor Preparation
The template for the M2 model construction was the high-resolution NMR structure of M2 channel protein (PDB code of 2RLF) [14]. In the M2 model construction of the 2009-H1N1 virus, the M2 protein sequence was taken from Genbank (accession number GQ385303), isolated in July 2009 from an H1N1 virus strain in Japan [10,15] (cf. Figure 5). Most synthetic inhibitors of M2 channel proteins, amantadine scaffolds based on adamantane derivatives were selected from published work by Gayday et al. [9], Du et al. [10], Eleftheratos et al. [18], Papanastasiou et al. [21], Tanner et al. [22], de Clercq [23], Tataridis et al. [24], Stamatiou et al. [25], Balannik et al. [26] and the others were newly created. These ligands have not been investigated for docking with M2 channel proteins either for the inside or outside positions before. The design principle of potential target drug-resistant influenza A on M2 mutants aimed to introduce an additional functional group on the aamantadine scaffold, including adding more amine, hydroxyl, linear, cyclohexane, aromatic and benzene groups, in order to increase the binding affinity of M2 inhibitors. The amino group in amantadine is likely the pharmacophore and is necessary to block hydrogen ion transport. Consequently, the scaffold was based on the adamantyl group and selection of suitable functional groups attaching to the scaffold to form new compounds that will potentially bind well to M2 channel proteins at the key residues in the binding site [27]. Prior to molecular docking simulation, all 200 new potential M2 inhibitors were built using Gaussview 4.1 software [28]. All compounds were optimized geometrically using Discovery Studio 3.0 Visualizer [29].

Virtual Screening with Autodock 3.05
In the molecular docking simulations, two types of M2 channel proteins representing H3N2 and 2009-H1N1 virus strains were used as receptor with two different docking positions inside and outside the channel. AutoDock Tools version 1.5.4 [30] was used to add polar hydrogens, assign Kollman charges and create grid binding boxes. The volume of each grid box was 30 × 30 × 30, with the default 0.375 Å spacing. For each type of M2 channel protein structures, AutoGrid version 3.05 was used to create affinity grids which centered at the active sites. Based on the coexistence of two possible binding positions and two existing M2 channel proteins [20], docking inhibitors into M2 channel will be separated into four categories: (1) docking inside the H3N2 M2 channel protein; (2) docking inside the 2009-H1N1 M2 channel protein; (3) docking outside the H3N2 M2 channel protein; and (4) docking outside the 2009-H1N1 M2 channel protein. The positive control for docking was obtained by re-docking amantadine and rimantadine extracted from the NMR structure [14] to locate the binding site on M2 channel protein before docking for 200 ligands. The binding box inside M2 channel proteins was positioned to encompass all three possible binding sites, namely Ala30, Ser31 and Gly34 and the binding box outside M2 channel proteins was set up to cover two important functional residues: channel gate Trp41 and channel lock Asp44 [15]. The binding affinities were calculated by AutoGrid version 3.05 using the following atom types: aromatic carbon (A), carbon (C), fluorine (F), iodine (I), nitrogen (N), hydrogen bond accepting N (NA), oxygen (O), hydrogen bond accepting O (OA), sulfur (S), hydrogen bond accepting (SA), silicon (Si), hydrogen (H) and electrostatic (e) [30]. The ligand set includes 200 new M2 inhibitors created from previously published studies and the others were newly created by modifying or adding more functional groups to the adamantane scaffold. AutoDockTools 1.5.4 was also used to merge nonpolar hydrogens, add charges and visually set up rotatable bonds for each ligand via AutoTors.
A Lamarckian genetic algorithm [31] was used to perform the docking experiments on AutoDock 3.05. The parameters were optimized as follows: trials of 100 dockings, population size of 50, random starting position and conformation, translation step range of 2.0 Å, rotation step range of 50 degrees, maximum number of generations of 27,000, elitism of 1, mutation rate of 2%, crossover rate of 80%, local search rate of 6%, 100,000 energy evaluations and docked conformations were clustered with the tolerance of 1.0 Å RMSD. Docking results were sorted by the lowest binding energy of the most populated cluster using AutoDockTools version 1.5.4. The top 10 hits from each of group binding inside and outside both H3N2 and 2009-H1N1 M2 channel proteins were selected for further analysis.

Pharmacophore Modeling
The top ten bound compounds were selected based on the lowest free energy of binding for pharmacophore analysis to give important insight into interactions between the top hits among M2 channel-inhibitors and common pharmacophore features in each group binding inside and outside. For this purpose, structure-and ligand-based pharmacophore modeling studies were carried out using the LigandScout 3.01 [32] and HipHop module of Discovery Studio 2.5 software [33], respectively. LigandScout generates the structure-based pharmacophore model based on the relevant interactions between the protein-ligand whereas Hip-Hop mainly focused on the possible common features present in the set of inhibitors.
3.3.1. Generation of Structure-Based Pharmacophore Models Using LigandScout 3.01 The top ten compounds binding on the inside and outside the M2 channel proteins with the lowest binding energy were used to generate the structure-based pharmacophore models [34]. The M2 channel-inhibitor observations were verified to compare the interactions between binding inside and outside of M2 channel proteins. The ligand interactions with critical amino acids present in the active site of M2 channel proteins pharmacophore based on best result of virtual screening provide a sufficient input to generate the structure-based. LigandScout was used to study the interactions between the M2 inhibitors and the amino acids in the two binding sites of M2 channel, as well as a tool for automatic construction and visualization of structure based pharmacophore model. LigandScout extracts and interprets ligand-receptor interactions such as hydrogen bond, charge transfer, hydrophobic regions of their macromolecular environment. Chemical features include hydrogen-bond donors and acceptors as directed vectors, positive and negative ionizable areas as well as lipophilic areas represented by spheres. In order to increase selectivity, excluded volume spheres are added to reflect potential steric restrictions. The 3D coordination of the interaction was obtained and resulted in specific interaction model that are able to map the ligands in their bioactive conformation. As a result, from the top 10 compounds binding at both sides, the most important inside interactions that can hold and stabilize the drug inside the M2 channel proteins were selected and visualised. The identification of important common chemical features from the top binding compounds inside and outside M2 channel proteins should be helpful to discover potent candidates to inhibit both the H3N2 and 2009-H1N1 virus strains. The significance of pharmacophore models mostly depends on the quality of the molecule structures used in generation of the pharmacophore conformation [35]. In this work, the training set molecule was selected from two different groups: the top 10 compounds binding inside the M2 channel proteins and the top 10 compounds binding outside the M2 channel proteins. The bond orders of these inhibitors were checked and verified before the generation of a pharmacophore model. The conformational flexibility of selected ligands was accomodated by creating multiple conformers to cover all representatives over a specified energy threshold. The HipHop module of Discovery Studio 2.5 software was employed to build a plausible binding hypothesis and identify a set of chemical features common to the most potent training molecule structures [36]. The general chemical features that considered in this training set were hydrogen bond donors and acceptors (HBDs and HBAs), and aliphatic hydrophobes. Eighteen hypotheses were generated in the outside group and seventeen hypotheses were created in the inside group, and the best hypothesis from each group was selected based on the high number of compounds fitting it as well as the high fit value of the hypothesis.

Conclusions
In this study the top highly potent anti-influenza A compounds from among 200 new compounds based on the amantadine scaffold were identified by a combination of different computational methodologies. The results also confirmed that the binding inside position was more favourable and stable than binding outside the M2 channel proteins. The results obtained based on virtual screening revealed that the top 10 compounds binding at both positions inside and outside the channel have higher binding affinity to both H3N2 and 2009-H1N1 M2 channel proteins than amantadine and rimantadine. Detailed pharmacophore analysis for the top hits among the M2 channel inhibitors also revealed the nature of interactions between functional groups of the top binding candidates with the M2 channel proteins. From this study, two compounds, I5 and I9, were predicted to be the most active inhibitors against both influenza A virus strains H3N2 and 2009-H1N1. However, due to the limitations of the method employed, such findings and additional mutation of M2 channel protein studies would require more accurate calculations to be confirmed. Molecular dynamics simulation should be done in the future to capture the flexible interactions between ligand-M2 channel proteins, together with in vitro and in vivo studies on the top binding compounds.