In Silico modeling of voltage-gated sodium channel alpha subunit to understand insecticide binding simulation in mosquitoes

The Voltage Gated Sodium Channel (VGSC) is critical for binding of different insecticides and play key role in insecticide resistance. An important mechanism of resistance to DDT and pyrethroids is termed knockdown resistance (kdr), caused by mutations in IIS6 domain of sodium channels. To attain a better management strategy for insecticide resistance and screening of new insecticide molecules, it is important to understand the three-dimensional structure of insecticide-binding domain of VGSC and its molecular interaction with insecticides. We constructed a theoretical model of ion transport domain–II of VGSC from mosquitoes, Culex quinquefasciatus. The stereochemistry of the model shows 91.1% residues are in the most favored region. Docking studies with DDT and deltamethrin indicated that deltamethrin showed interaction with Thr929, Met918, Ile936, Cys933, Leu925, Glu881, Met857 and Gly866 and DDT showed interaction with Ile936, Thr929, Ser878, Phe863, Gln864, Trp861 and Met857. We also predicted that mutation of Thr929 should confer resistance to both DDT and deltamethrin.


Introduction
The Voltage Gated Sodium Channel (VGSC) mediates the initiation and propagation of action potential in the nervous system and other excitable cells. The critical role of the sodium channel in membrane excitability has made it the target for a variety of neurotoxins during evolution (Wang and Wang, 2003). During an action potential, the sodium channel undergoes transactions between closed-resting, activated, and inactivated functional states and neurotoxins alter various channel properties, including ion conductance, ion selectivity, activation, or inactivation. "There are at least ten separate binding sites for ligands on the sodium channel, including those for local anesthetics and anti-convulsants" (O'Reilly et al., 2006). The site classified as '7' has been recognized as binding site for the DDT [1,1,1-trichloro-2,2-bis (p-chlorophenyl) ethane] and pyrethroid insecticides (Wang and Wang, 2003).
Prior to the 1970s, DDT was used intensively in agriculture and in vector control, but following concerns over its environmental impact, its use in agriculture was discontinued or banned, but it is still recommended in vector control programs by World Health Organization (WHO). Pyrethroid, the synthetic analogues of the naturally occurring pyrethrum from the flower extracts of Chrysanthemum species, is a favored alternative of DDT. Because of the relatively low mammalian toxicity, low persistence in the environment and high excito-repellency activity, pyrethroids represent a major class of insecticide used to control many agriculturally and medically important arthropod pests.
Lymphatic or bancroftian filariasis is considered as the predominant infection in the continental Asia (Gyapong et al., 2005). The mosquitoes, Culex quinquefasciatus is the principal vector of the parasitic worm Wuchereria bancrofti the agent of bancroftian filariasis throughout the continental Asia. Due to indiscriminate and intensive use of insecticides, many insects including mosquitoes have developed resistance to these compounds. An important mechanism of resistance to DDT and pyrethroids is termed as knockdown resistance (kdr) was first discovered in houseflies (Milani, 1954) and subsequently in many other insect and arachnid species (Soderlund and Bloomquist, 1990;Soderlund and Knipple, 2003). Extensive research has shown that kdr or kdr-like mechanisms are resulted by mutations in sodium channels (Dong, 2002;Soderlund and Knipple, 2003;Vais et al., 2001). Like mammalian sodium channel alpha subunits, the primary structure of insect sodium channel proteins consists of four homologous domains, each containing six transmembrane segments (S1 -S6) (Loughney et al., 1989). Mutations in the domain-II region of the channel protein are commonly responsible for insecticide resistance. The most common mutation is leucine to phenylalanine at IIS6 domain and is referred to as the kdr mutation.
To attain a better management strategy for insecticide resistance, development of effective interference mechanisms and screening of new insecticide molecules, it is essential to understand the three-dimensional structure of voltage-gated sodium channel and the molecular interaction between the target site and its ligands. Therefore, we designed a three-dimensional theoretical model of the insecticide-binding domain (i.e., IIS6 domain) of voltage-gated sodium channel alpha subunit from Culex quinquefasciatus and studied the molecular interactions between channel protein and insecticides. We computed this molecular interaction using two insecticides, DDT, and deltamethrin. Automated docking studies were used to perform this insecticide binding simulation analysis.

In Silico modeling
The amino acid sequence of Voltage gated sodium channel alpha subunit of Culex quinquefasciatus (Taxonomic Identifier: 7176, NCBI) was retrieved from the Uniprot database (http://www.uniprot.org/uniprot/A5I9E7; accession number: A519E7). In this sequence, the site of kdr mutation that occurred due to substitution of leucine to phenylalanine is at position 1016, whereas in general this mutation found at position 1014 in other referral sequences (Martinez-Torres et al., 1998, 1999Chandre et al., 1998;Wondji et al., 2008). However, we did not make any manual adjustment to the retrieved sequence during modeling. But during discussion of ligand binding simulation, this error had been taken into consideration and based on the available sequence in literatures, the position of all amino acid residues were accepted as -2 (minus two) from their original position in the three-dimensional model.
The sequence was submitted to the Pfam to search for its families and domains (http://pfam.sanger.ac.uk/search). BLAST search was used to find the homology of the query sequence with other known sequences of the database. Initially full-length protein sequence (accession number: A519E7) was submitted to the server 3D-JIGSAW (version 3.0) POPULUS (http://bmm.cancerresearchuk.org/~populus/populus_s ubmit.html) in search of suitable templates for homology modeling of VGSC protein. The server identified templates using HMM (Soding, 2005) and returned alignments were used to build the model. All models were preselected using POPULUS ENERGY, gaps and missing residues were closed and filled using POPULUS REPAIR, finally all models were recombined using basic POPULUS approach (Offman et al., 2006), where Genetic Algorithm (GA) was used for conformational space search engine. All models were ranked using a fine and coarse energy function weighted according to the highest sequence identity. The server constructed several models on split site basis using different templates rather than building a single model covering the whole sequence.
After analyzing all the initial models, we did not found any suitable template for construction of the whole VGSC protein model. Therefore, depending on the sequence id and energy function, only domain-II of VGSC (residues 833 to 1047) was modeled. We selected the X-ray crystal structure of the rat brain Kv 1.2-Kv 2.1 Paddle Chimera Channel (PDB Code: 2R9R, chain 'B'), a member of the Shaker family of voltage depended potassium channel as a structural template to generate a homology model of VGSC domain-II. The model was then recombined using basic POPULUS to other 29 models (Table 1), which were selected according to their coverage, sequence id and domains, for the construction of the final three-dimensional structure of VGSC domain-II (residue coverage: 833 to 1047). Note: The 3D-JIGSAW POPULUS server while used in interactive mode generates above table The backbone conformation and stereo-chemical properties of the constructed model was evaluated by Psi/Phi Ramachandran plot (Ramachandran et al., 1963) obtained from PROCHECK (Laskowski et al., 1993). The constructed model of VGSC was submitted to NIH MBI Laboratory for Structural Genomics and Proteomics Sever, an online web interface (http://nihserver.mbi.ucla.edu/) which provides validation report from PROCHECK. STRIDE (http://webclu.bio.wzw.tum.de/cgi-bin/stride/stridecgi.py) which uses hydrogen bond energy and main chain dihedral angles to recognize helix, coil, and strands was used to predict the secondary structure of the modeled protein. The model was submitted to ConSurf (http://consurf.tau.ac.il/) an automated web based server for the identification of functional region in proteins of known three-dimensional structures by estimating the degree of conservation of the amino acid, based on phylogenetic relationship between their close sequence homologues. Given the 3D-structure of the protein as an input, the ConSurf server extracts the sequence from the PDB file, automatically performs a search for close homologous sequences of the protein of known structure using PSI-BLAST, and then aligns them. The multiple sequence alignment is then used to build a phylogenetic tree. Then the program calculates the conservation scores using either an empirical Journal of Mosquito Research 2015, Vol.5, No.23, 1-8 http://jmr.biopublisher.ca Bayesian (Mayrose et al., 2004) or the Maximum Likelihood method.
The three dimensional structure was also submitted to Meta-Function Signature (MFS) server (http://protinfo.compbio.washington.edu/mfs/) for quantitative measures of functional importance of each constituent amino acids, calculated by combining sequence, structure, evolution and amino acid property information (Wang et al., 2008). Presence of pockets in the modeled protein was predicted using Computed Atlas of Surface Topography of Proteins (CASTp) server (http://sts-fw.bioengr.uic.edu/castp/).

Ligand docking
The chemical structure of deltamethrin and analogues of DDT, having more than 90% structural similarity with the DDT, were retrieved from the PubChem Database (http://pubchem.ncbi.nim.nih.gov/). Following the reference of Rajesh et al. (2007), the 1-chloro-4-[1,2,2,2-tetrachloro-1-(4-chlorophenyl)-eth yl was selected as the best ligand for docking. The structure of deltamethrin and DDT analogue was retrieved into two-dimensional Structure Data File (SDF) format. The three-dimensional coordinates in PDB format were generated using the online SMILES Translator and Structure File Generator (http://cactus.nci.nih.gov/services/translate/) program. Docking was performed using Patch Dock Server, an online web server (http://bioinfo3d.cs.tau.ac.il/PatchDock/index.html). Default parameters of Patch Dock program were used for docking purpose. The server returned 500 docking predictions, from which 10 top solutions, selected by the server were retrieved. From those 10 solutions, three docking results were chosen that had binding contacts with residues known to be important in insecticide resistance and where the ligand came within 4Å of these residues. Figures (2, 7 and 8) were produced using the PyMOL molecular graphics system (DeLano Scientific, San Carlos, CA, U.S.A.).

Modeling of domain-II of voltage gated sodium channel
The protein sequence (accession number: A519E7) retrieved for this study was 2129 amino acid in length with molecular weight of 238,442 Da. Pfam search result showed 20 Pfam-A matches to the sequence (6 significant and 14 insignificant). The Pfam results further reveled that there were five domains in the query sequence, four of which belongs to ion transport domains (PF00520) and one sodium transport associated channel (PF06512).
We performed the PSI-BLAST that failed to find a suitable template for modeling the whole sequence. Therefore, only domain-II of VGSC, important for insecticide binding, was modeled. Figure 1 shows the amino acid sequences highlighted with grey color used for construction of the model. The leucine residue (red in color) highlighted with yellow color is the site for kdr mutation (Fig. 1). Table 1 presents the templates that were used for POPULUS recombination to generate the final VGSC domain-II model. The server returned top five models for the query sequence and we selected the final model depending on the sequence id and energy (-249.96 kcal/mol) function (Fig. 2). Figure 1 The amino acid sequence of voltage-gated sodium channel alpha subunit of Culex quinquefasciatus, retrieved from the Uniprot database (accession number A519E7). Sequence coverage and kdr mutation point in the model are shown in grey and yellow color respectively PROCHECK estimated the stereo-chemical quality of the modeled VGSC protein (Fig. 3). The phi/psi angles of 91.1% residues (total 204 residues) fell in Journal of Mosquito Research 2015, Vol.5, No.23, 1-8 http://jmr.biopublisher.ca 5 most favored regions, 6.7% residues (total 15 residues) lied in the additional allowed regions, 1.8% residues (total 4 residues) fell in the generously allowed regions; only 0.4% residues (total 1 residue) of residue lied in the disallowed regions. Secondary structure assignment using STRIDE provided the physical features of the modeled structure, such alpha helix, 3-10 helix and coil / turn (Fig. 4). The ConSurf and MFS server was used to extract information about important residues, which are of functional value. The ConSurf provides evolutionary related conservation scores for residues. Conservation grades are projected onto the molecular surface of the protein model to reveal the patches of evolutionarily highly conserved residues that are often important for biological function (Fig. 5). Quantitative measures of functional importance of each constituent amino acid was calculated by combining sequence, structure, evolution and amino acid property information using Meta-Function Signature (MFS) program. MFS server predicted top 10 highest scoring residues are Glu851 (0.86), Lys855 (0.88), Tyr864 (0.91), Asp873 (0.97), Ser880 (0.82), Arg899 (0.92), Arg902 (0.97), Lys905 (0.82), Glu987 (1.00), Trp988 (0.95). CASTp program illustrated the presence of 32 pockets in the modeled protein, which are important for ligand interaction (Fig. 6). Figure 3 Ramachandran plot of our modeled VGSC domain -II generated by PROCHECK program. Based on an analysis of 118 structures of resolution of at least 2.0 Angstroms and R-factor no greater than 20%, a good quality model would be expected to have over 90% in the most favored regions. The plot statistics are: residues in most favored region [A, B, L] -204 (91.1%); residues in additional allowed regions [a,b,l,p] -15 (6.7%); residues in generously allowed regions [~a,~b,~l,~p] -4 (1.8%); residues in disallowed regions -1 (0.4%); number of non-glycine and non-proline residues -224 (100.0%); number of end-residues (excl. Gly and Pro) -1; number of glycine residues (shown as triangles) -14; number of proline residues -8; Total number of residues -247 The hypothetical docking studies on the theoretical model of domain-II of VGSC with two types of insecticides was performed to explore the importance of voltage gated sodium channel residues involved in the interaction with these insecticide. In Silico docking studies were performed with DDT, an organochlorine and deltamethrin, a synthetic pyrethroid. Both DDT and pyrethroids exert a similar functional effect on the insect voltage gated sodium channel, namely stabilization of the open channel state, slowing of repolarization after an action potential, and repetitive discharges (O'Reilly et al., 2006). In addition, kdr mutation may confer resistance for both insecticides simultaneously (cross-resistance). These observations suggest that DDT and pyrethroids may share an overlapping binding site. Docking results indicated that deltamethrin showed atomic interaction with Thr931, Met920, Ile938, Cys935, Leu927, Glu883, Met859 and Gly868 (Fig. 7). On the other hand DDT showed interaction with Ile938, Thr931, Ser880, Phe865, Gln866, Trp863and Met859 (Fig. 8). The amino acids of VGSC domain-II are colored in the range from turquoise to maroon based on conservation grades according to ConSurf description. The VGSC-IIS6 protein is represented as a spacefill model, where the residue conservation scores are color-coded onto its Van der Waals surface. The color-coding bar shows the coloring scheme; conserved amino acids are colored maroon, residues of average conservation are white, and variable amino acids are turquoise As mentioned in the methodology section, the site of kdr mutation in the sequence, retrieved for the present study is at position 1016, whereas in general this mutation found at position 1014 in other referral sequences. Thus, the positions of all amino acid residues in the three-dimensional model were accepted as -2 (minus two), based on other available sequences in the literatures (i.e. Thr931, Met920, Ile938, Cys935, Leu927 should read as Thr929, Met918, Ile936, Cys933, Leu925 respectively). In this connection, our docking predictions support the findings of Usherwood et al. (2007) who found  substitutions M918T, L925I, T929I, C933A decreased   deltamethrin potency, M918T, L925I, T929I decreased  permethrin potency and T929I, L925I, and I936V decreased fenfluthrin potency in Drosophila. Our docking prediction also support those of Usherwood et al. (2007) in the conclusion that DDT potency could badly affected by T929I and reduced by L925I, L932F and I936V substitutions in Drosophila. In our docking study, we found that both DDT and deltamenthrin interact with Thr929. Therefore, we also predicted that any mutation of Thr929 should confer resistance to both DDT and deltamethrin in mosquitoes. Previous studies on this mutation where a substitution of Thr to Ile in diamond-back moth (Plutella xylostella) have shown that this is indeed the case, DDT being completely ineffective against larvae of the pyrethroid resistant Plutella strain (Schuler et al., 1998). Figure 6 Available pockets for ligand interaction in the VGSC domain-II as predicted by CASTp program. For better visualization, annotated pockets are presented with corresponding sequence map, where residues in highlighted pocket are highlighted in the same color as in the structural visualization However, the most common resistance mutation L1016F, which gives the kdr phenotype, does not however appear to be a binding determinant, as it is located far apart from the main binding pocket in the model. This result is consistent with the finding of O' Reilly et al. (2006) and Rajesh et al. (2007). However, the mutation at this point (L1016F) lower the affinity of open channels for pyrethroids by 10-30 fold, and decrease the availability of open channel states due to enhanced closed-state inactivation, thereby limiting the number of high affinity binding sites available for pyrethroids (Vais et al., 2000). However, the docking results of this study should be interpret with caution due to the fact that docking was done on an isolated domain II model not the whole pore region of the channel. A whole VGSC model of mosquito is needed for better understanding of these molecular interactions. The current study provided a three dimensional structure of domain-II of voltage-gated sodium channel alpha subunit in Culex quinquefasciatus which is very essential for binding of insecticides. The results of insecticide binding simulation may be helpful in future experimental research on site directed mutagenesis for insecticide resistance study. The results of insecticide binding simulation also identified the possible new targets for investigation, which could help in explaining the insecticide resistance mechanisms against DDT and pyrethroid.

Ligand binding simulation
Hence, the current investigation is very important at both fundamental and applied level.