Malaria is one of the major health problems of humans caused by parasites belonging to the genus Plasmodium.¹ According to WHO, there were an estimated 241 million malaria cases worldwide and 627000 deaths in 2020.² In the South-East Asia Region, three countries accounted for 98% of the total reported cases, the main contributor being India (58%), followed by Indonesia (30%) and Myanmar (10%). According to the World Malaria Report, India represents 3% of the global malaria burden.³ Use of insecticides such as pyrethroids (e.g., permethrin, deltamethrin and lambda-cyhalothrin), organochlorines (e.g., DDT), carbamates (e.g., bendiocarb, propoxur) and organophosphates (e.g., malathion, fenitrothion, pirimiphos-methyl) is the primary strategy for controlling mosquito vectors.^4,5 However, there are reports of insecticide resistance developed by vectors.⁶ Studies have observed that mutations in the Voltage-Gated Sodium Channel (VGSC) are associated with insecticide resistance, termed knockdown resistance (kdr), in Anopheles mosquitoes.¹ Sodium channels are the primary targets of several classes of synthetic insecticides, including DDT and permethrin. The binding of permethrin causes prolonged polarization of the channel protein, stabilizing open-state leading to paralysis and death of the mosquitoes. Mutation in kdr gene is found to be the major reason of permethrin insecticide resistance in mosquito vectors.^7–9 In neurons and other excitable cells, voltage-gated sodium channel originates and disseminates action potentials. They propagate along the muscle fibres, somatodendritic compartments of neurons and axons of nerves.¹⁰

A wide range of applications in life sciences research is influenced by the 3D structures of proteins, which offer valuable insights into their molecular function.¹¹ The degree to which a protein satisfies the requirements for both essential amino acids and physiological needs can be studied by the overall quality of the protein.¹² The docking methodology is a standard computational tool in drug design for lead compound optimization and virtual screening to find novel bioactive molecules. It can also predict the experimental binding modes and affinities within the binding site of specific receptor targets.¹³ According to the overall energy of the system, docking simulations can forecast an optimum docked conformer.¹⁴ Several studies have reported the mutational changes in the kdr gene and the development of insecticide resistance in Anopheles mosquitoes. However, whether the kdr gene mutation and structural changes in VGSC affect insecticide binding affinity has not been explored. Therefore, the present study investigated the structural alteration of VGSC protein and its binding affinity to permethrin insecticide.

The amino acid sequences of VGSC of wild-type Anopheles mosquito (Anopheles funestus) (Accession no. AUA17903) and mutated sequence (L1014C)¹ were retrieved from NCBI¹⁵ (https:// www.ncbi.nlm.nih.gov). Both the sequences were submitted to SWISS-MODEL (https://swissmodel.expasy.org) to build the 3D structures.¹⁶ All the models were downloaded in PDB format.

The overall structure of protein model was analyzed by Procheck,¹⁷ QMEAN (Qualitative Model Energy Analysis),¹⁸ ANOLEA (Atomic Non-Local Environment Assessment),¹⁹ ERRAT,²⁰ and Ramachandran Plot (https://swissmodel.expasy.org). The overall quality of the models was studied in UCLA-DOE LAB-SAVES v6.0 (http://saves.mbi.ucla.edu). The FASTA sequence was submitted to InterPro²¹ to analyze the protein by classifying it into families and domains. Protein Secondary Structure was predicted by Phyre2.²² Both the wild-type and mutant proteins were added to Chimera Software²³ to study the structural similarities. The RMSD value was calculated by putting the proteins in PyMol software.

The 3D structure of wild and mutant VGSC proteins was converted into .pdbqt file using OpenBabel.²⁴ The 3D conformer of Permethrin (CID: 40326) was acquired from the PubChem database (http://pubchem.ncbi.nlm.nih.gov/) in the SDF file format and converted to .pdbqt format using the software OpenBabel.²⁴ The active pocket of the model proteins was created based on the ligand interaction data of the template protein. The grid parameters were set as 48, 64, 92 for axes X, Y and Z and 166.006, 181.239 and 199.338 as center-coordinates, respectively. Molecular docking was allowed keeping the default exhaustiveness and energy range at 8 and 4, respectively. The docking was carried out in triplicate (n=3) using Autodock Tools-1.5.6²⁵ and the output was visualized in Discovery Studio-2020.

The amino acid sequence of VGSC contains 2116 residues. After submitting the amino acid sequences in SWISS MODEL server, the system generated four numbers of models. The best model was selected based on QMEAN value (Figure 1). Using the QMEAN tool, two different QMEAN Z-score values were determined. The QMEAN value for wild-type protein is -6.05 and the value for mutated protein is -6.09 (Figure 2). Figure 1d showed the superimposition of both the wild-type and mutated VGSC proteins. Template for wild type protein P35500.1.A Sodium channel protein para AlphaFold DB model of SCNA_ DROME (gene: para, organism: Drosophila melanogaster (Fruit fly). The root means square deviation (RMSD) of wild-type and mutated protein backbone was found to be 0.117Å.

Figure 1:
3D structures of (a) template protein (b) wild-type protein (c) mutated protein and (d) superimposition between wild and mutated ion channel proteins (wild protein-cyan and mutant-yellow).

Figure 2:
Protein quality analysis in Structure Assessment tool of SWISS-MODEL (a) Local Quality Estimation of wild-type protein (b) Local Quality Estimation of mutated protein (c) Normalized QMEAN of wild-type protein (d) Normalized QMEAN of mutated protein (e) Coloured sequence of wild-type protein (f) Coloured sequence of mutated protein.

The protein conformation of both the wild and mutant VGSC was studied using a Ramachandran plot. Figure 3 shows an almost similar distribution of amino acid residues in the Ramachandran plot. The core region of wild-type and mutated protein has 85.57% and 84.95% of distribution, respectively. The outliers of wild-type and mutated protein show 4.53% and 4.31% residues distribution. The rotamer outliers of wild-type and mutant protein have been found to be 1.72% and 2.09%, respectively. The overall quality value of wild-type and mutated VGSC protein was found to be 85.204 and 87.541, respectively. The overall G-factor of the wild-type protein is -0.18 and the Mutated protein is -0.21 (Figure 4).

Figure 3:
Ramachandran Plots of (a) wild-type protein and (b) mutated protein.

Figure 4:
Overall Quality factor of (a) wild-type protein and (b) mutated protein obtained from ERRAT program of SAVES v6 server.

The domain studies have revealed that the protein contains five numbers of domains. The amino acid residue numbers 151 to 437 and 1857 to 1892 form the EF-hand domain which has the molecular function of calcium ion binding. The residue numbers 542 to 723 form the cytoplasmic domain, which has the function of Na⁺ ion binding. The amino acid residues 793 to 1019, 1253 to 1412, 1416 to 1550 and 1600 to 1850 form the Ion transport domain. This domain has the molecular function of ion channel activity. The residue numbers 1029 to 1249 form the Sodium ion transport-associated domain which has the molecular function of voltage-gated sodium channel activity. The amino acid residue numbers 1541 to 1593 form the VGSC-alpha subunit, inactivation gate domain. Most electrically excitable cells use this region to produce the depolarizing phase of the action potential. They can exist in three different states: resting, in which the channel is closed; activated, in which the channel is open; and inactivated, in which the channel is closed and resistant to opening.

The secondary structures of the domain were obtained from the Phyre2 server. Domain 1 is formed by two groups of sequences which contain 287 and 37 residues where 68% of Alpha helix, 0% of Beta strands, 51% of TM helix and 5% of disordered residues were observed. Domain 2 forms by two groups of sequences which contain 43 and 182 residues where 7% of Alpha helix, 23% of Beta strands and 51% of disordered residues were observed. Domain 3 forms by three group of sequences which contains 222, 160, 136 and 251 residues where 82% alpha helix, 0% beta strand, 64% tm helix and 3% Disordered residues were observed. Domain 4 contains 66 residues where 58% of alpha helix, 4% of beta strand and 47% of disordered residues were observed. Similarly, Domain 5 contains 52 residues where 58% of alpha helix, 4% of beta strand and 40% of disordered residues were observed.

To further understand how the ion channel resistance occurs to permethrin insecticide, molecular docking was conducted. Figure 5 shows the docking of the wild-type and mutated VGSC and their 2D display indicating the amino acid involved in the interactions. The binding affinity of permethrin was found to be stronger with wild-type protein (docking score -7.9 kcal/mol). A total of 16 amino acid residues (GLN385, CYS973, GLY974, ASN1003, LEU1007, ALA1486, THR1487, PHE1488, LYS1489, SER1528, PHE1529, LEU1532, ASN1533, THR1780, PHE1828 and LEU1829) were found involved between wild-type protein and permethrin. Out of 16 amino acids, GLY974, ASN1003, and THR1780 showed H-bond interaction, LEU1532 showed Pi-sigma interaction. The amino acid residue LEU1007 and PHE1828 form the Alkyl interaction. Remaining residues GLN385, CYS973, ALA1486, THR1487, PHE1488, LYS1489, SER1528, PHE1529, ASN1533, and LEU1829 involves in Van der Waals interaction. Mutated VGSC protein showed slightly weaker binding affinity with permethrin (docking score, -6.9 kcal/ mol). The interaction between mutant protein and permethrin showed that 15 amino acid residues, THR384, GLN385, SER415, PHE416, VAL419, ASN420, LEU423, GLY974, ASN1003, LEU1004, ASN1008, ILE1536, PHE1828, ILE1832 and ILE1836 were involved in the interaction. Residues ASN420 and ASN1003 are involved in H-bond interaction, VAL419 and PHE1828 in Pi-sigma and Pi-alkyl interactions, respectively. The remaining residues, THR384, GLN385, SER415, PHE416, LEU423, GLY974, LEU1004, ASN1008, ILE1536, ILE1832 and ILE1836 were involved in Van der Waals interactions.

Figure 5:
Structure showing (a) 3D Complex of wild-type ion channel protein and Permethrin (b) Interaction of wild-type protein and Permethrin (c) 3D Complex of mutated ion channel protein and Permethrin (d) Interaction of mutated ion channel protein and permethrin obtained from docking, visualized in Discovery Studio.

Ion channels are integral membrane proteins that regulate the movement of ions across the cell membrane. They are involved in several biological processes and malfunction of which can lead to many diseases as well as the death of the organism. The effectiveness of malaria control programs is compromised by insecticide resistance in Anopheles mosquitoes.²⁶ The increasing resistance ability of mutated VGSC resulted in resistance of permethrin, a commonly used insecticide of malaria-causing vectors, by inhibiting their function. By reviewing the literature worldwide, it has been observed that very few studies were done on insecticide resistance of the VGSC of Anopheles mosquitoes. This study built the 3D structure of both wild-type and mutated voltage-gated ion channels. After observing the qualitative model energy analysis, it has been observed that the mutation at 1014 codon changes the conformation of VGSC protein leading to the decrease in the overall QMEAN of the mutant protein. Good agreement between the model and experimental structures of comparable size is indicated by QMEAN z-scores near zero.¹⁸ The Ramachandran plot has served as the cornerstone of protein structural validation. With the availability of more and more experimentally determined models of protein 3D structures, their intricate structure has been continuously analyzed and improved.¹⁷ The present study observed a slight difference in terms of the distribution of amino acid residues in the Ramachandran plot. The mutated protein showed a slightly higher percentage of residues in the disallowed region of the plot. Similarly, the wild-type VGSC protein of Anopheles gambiae showed 88.6% residues in the favoured region.²⁷ Similarly, 93.9% of the residues fell in the most favoured regions in PR-1a Protein of Hordeum vulgare.²⁸ In this study, the Ramachandran favorued residues of wild-type protein were 85.57% residues; in mutated protein, the Ramachandran favoured was 84.95%. The more favoured percentage residues indicate a more stable protein structure.

In this study, the overall G-factor of wild-type protein is -0.18 and the mutated protein is -0.21. The dihedral angles and main-chain covalent forces had an average G-factor of -0.21 in PR-1a Protein of Hordeum vulgare, which was higher than the threshold value of -0.5.²⁸ The overall determining quality factor of ERRAT was 94.40 in PR-1a Protein of Hordeum vulgare.²⁸ In the VGSC of Anopheles, the ERRAT value of wild-type is 85.204 and mutated is 87.541, stating that the error values of each residue are quite small. Molecular docking studies showed a stronger binding affinity of permethrin with wild-type channel protein compared to mutated protein. It acts on the nerve cell membrane and disrupts the sodium channel leading to the depolarization of the membrane and finally death of the mosquitoes.^29,30 Several studies have reported the importance of ion channel as major drug target. Philanthotoxins from wasp venom provide useful pharmacological tools acting on insect and vertebrate ligand-gated ion channels. Kachel et al.³¹ reported strong binding affinity of philanthotoxins isolated from wasp against vertebrate ligand-gated ion channels. Recent studies by Qiao et al.³² reported the ion-channel inhibitory property of Azobenzene-isoxazoline. Similarly, isocycloseram, a novel insecticide has also been found to possess strong inhibitory activity against GABA-associated ion-channel.³³ In the present study, we observed slightly weaker binding affinity of permethrin insecticide with mutated ion-channel protein compared to wild non-mutated protein. The present study therefore suggests that the mutations in ion channel leads to the structural alterations of the protein thereby decreasing the binding affinity of permethrin.

Anopheles is the core vector of malaria disease. Insecticides that suppress the malaria mosquito are essential for controlling the disease. Due to the intense selection pressure that the pyrethroids are applying, changes have occurred in populations of the same species and populations of distinct species and they are dispersing swiftly. However, according to recent findings, the VGSC is quickly acquiring point mutation resistance to insecticides. Using the in silico molecular modeling and docking simulation, the interactions of permethrin with the VGSC proteins from Anopheles that were both unaltered and altered were detected. The study of the molecular docking data revealed that the amino acid residues of VGSC change the molecular interaction after mutation. This inference suggests that mutation in codon 1014 affects the complex stability, effectively resisting the effect of the insecticide permethrin. This investigation could be utilized to understand better the mechanism of insecticide resistance property of transmembrane protein against different commercially available insecticides. However, a molecular simulation study can be done to know the detailed insecticide resistance mechanism in Anopheles.

The authors would like to acknowledge the facilities provided by the Department of Zoology, Bodoland University for providing all the necessary facilities to carry out the work.