In silico identification and study of potential anti-mosquito juvenile hormone binding protein (MJHBP) compounds as candidates for dengue virus - Vector insecticides.

Dengue has become a huge global health burden. It is currently recognized as the most rapidly spreading mosquito-borne viral disease. Yet, there are currently no licensed vaccines or specific therapeutics to manage the virus, thus, scaling up vector control approaches is important in controlling this viral spread. This study aimed to identify and study in silico, potential anti-mosquito compounds targeting Juvenile hormone (JH) mediated pathways via the Mosquito Juvenile Hormone Binding Protein (MJHBP). The study was implemented using series of computational methods. The query compounds included pyrethroids and those derived from ZINC and ANPDB databases using a simple pharmacophore model in Molecular Operating Environment (MOE). Molecular docking of selected compounds' library was implemented in MOE. The resultant high-score compounds were further validated by molecular dynamics simulation via Maestro 12.3 module and the respective Prime/Molecular Mechanics Generalized Born Surface Area (Prime/MM-GBSA) binding energies computed. The study identified compounds-pyrethroids, natural and synthetic - with high docking energy scores (ranging from 10.91-12.34 kcal/mol). On further analysis of the high-ranking (in terms of docking scores) compounds using MD simulation, the compounds - Ekeberin D4, Maesanin, Silafluofen and ZINC16919139- revealed very low binding energies (-122.99, -72.91 -104.50 and,-74.94 kcal/mol respectively), fairly stable complex and interesting interaction with JH-binding site amino acid residues on MJHBP. Further studies can explore these compounds in vitro/in vivo in the search for more efficient mosquito vector control.

Introduction

Dengue has become a huge global health burden. It is currently recognized as the most rapidly spreading mosquito-borne viral disease. The disease is caused by any one of the four single- and positive-stranded dengue serotypes (DENV 1, 2, 3, 4) belonging to the genus and family of Flavivirus and Flaviviridae, respectively [[1], [2]]. The virus is wide spread throughout the tropics and presents with a spectrum of diseases which climaxes as severe dengue [3]. Transmission is primarily via Aedes aegypti and secondarily by other species within the Aedes genus, such as Aedes albopictus [1]. The prevalence and incidence of the virus have risen drastically in the last couple of decades. According to a modeling estimate in 2013, there was a startling statistics of 390 million global dengue infection per year, of which only 96 million manifested clinically and approximately 20,000 cases died [[2], [4]]. More recently, Zeng et al. [2] have also argued that the prevalence and incidence rate have been underestimated and thus, reported that from 1990 to 2017 (approximately three decades), the number of dengue cases rose from 23.3 million to 104.8 million (22.2%); and death cases from 16,975 to 40967 (40.46%). In Africa, the surveillance data have been quite poor, but strong evidences of past dengue virus exposure and infection are well documented in Nigeria [5], and outbreaks in Zanzibar, Burkina-faso, Egypt, Senegal and South Africa [6]. With this overwhelming statistics, it is imperative that intervention efforts and technique are scaled up. Unfortunately, till date, treatment of dengue is still limited to supportive care and symptomatic treatment. There is no anti-dengue drug or vaccine except for Dengvaxia approved by FDA for use in only serogative positive patients [2].

A typical good approach to controlling parasitic diseases is to target the vectors of the parasite – chemically or biologically. Quite a number of insecticides - implemented via indoor residual spraying and long lasting insecticide-treated nets - have been employed in vector control [7]. These include organochlorine, organophosphates, carbamates, pyrethroids and neonicotinoids [8]. Although, some of these insecticides have been relatively efficient, the use and success rate has been limited by the emergence of resistance to and reported toxicity for most of them [9]. The introduction of synthetic insecticides, and the emergence of insecticide resistance date back to 1940s, following the advent of a synthetic pesticide, dichlorodiphenyltrichloroethane (DDT). It is now evident that continual use of a class of insecticides or those with similar modes of action exert a selection pressure for resistant strains on vectors [[10], [11]]. So, it is important to identify and characterize more potent, less toxic insecticidal agents with a different mechanism of action.

Targeting or obstructing the normal function of the endocrine system of insects is recognized to be a logical and promising approach to the development of bio-rational vector control agents/management systems - huge successes have been achieved in this regard [12]. Insect metamorphosis is coordinately regulated by the action of two insect-specific hormones – (1) 20-hydroxyecdysone (20E), which triggers the successive molts throughout the life cycle, and (2) juvenile hormone (JH), which represses the transition to the adult stage [13]. Juvenile hormones (JHs) are sesquiterpenoids synthesized by the corpora allata (CA), and they play key roles in maintaining the larval stage of insects and transition to adult forms [14]. The signaling pathways mediated by JH encompass major pathways such as metamorphosis, cellular immunity, trehalose metabolism, reproduction, egg production, vitellogenesis [15]. Hence, JH signaling pathway may be considered a suitable target in the development of novel insecticides with low toxicity to non-target vertebrate organisms - they are specific to insects and other arthropods and do not exist in vertebrates. However, despite theses overwhelming pieces of evidences that posit JH endocrine system as a potential target of insecticides, there are only a few well known insecticidal agents (referred to as insect growth regulators (IGRs)) that target JH signaling pathways [16]. This may have been due to the fact that the JH target pathways, proteins, transport and synthesis accessories had not been adequately characterized. Many years after the discovery of JH by Sir Vincent B. Wigglesworth, there is still paucity of data on the molecular mechanism of this hormone and its associated proteins in specific arthropods, particularly mosquitoes [14]. The target molecules may include JH receptors/transcription factors (methoprene-tolerant protein; Krüppel-homolog 1 (Kr-h1)), JH binding proteins and enzymes involved in JH synthesis (Juvenile acid methyl transferase, JH esterase, epoxidase). The protein methoprene-tolerant (Met), which was first reported in Drosophila melanogaster (and conserved in insects and other arthropods) is recognized as the receptor for the JH ligand. It was identified as a member of the basic helix-loop-helix Per-ARNT-Sim (bHLH-PAS) family of transcription factors, which are essential for regulating critical gene regulators of gene expression networks underlying many essential physiological and developmental processes [13]. Krüppel-homolog 1 (Kr-h1) is an early inducible gene in the JH signaling pathway that activates the transcription of many JH-inducible genes. The JH synthetic pathway enzymes include: Juvenile hormone acid methyltransferase (JHAMT) - catalyses the transformation of farnesoic acid (FA) to methyl farneosate (MF); Epoxidase (epox) – catalyses the conversion of MF to JH; JH esterase (JHE) converts JH into JH-acid [17].

Juvenile hormone binding protein is a member of the odorant-binding protein (OBP) family and orthologs. They are present in the genomes of Aedes, Culex, and Anopheles mosquito species. Juvenile binding proteins specifically bind JH II and JH III but not eicosanoids and are thought to play a significant and essential role in the transportation of JH hormones within the insects [18]. They have been reported as targets for novel anti-mosquito insecticide [[19], [20], [21]]. A significant development in this area in recent times is the determination of the 3D x-ray crystallographic structure of this protein in complex with JH for mosquito species by [18].

Computer-aided molecular design is a rational approach that comprises a broad range of theoretical and computational approaches often employed in modern drug design [22,23]. It has become an essential tool for lead screening, optimization and design of new potent inhibitors, including insecticidal agents due to current advances of biochemistry and structural biology [24]. The rationale is basically to expedite the process and also reduce the cost relative to standard experimental approaches [25]. Additionally, they are amenable to working with huge, diverse databases of compounds. For a very long time, the detailed data on the mechanism of JH action, pathways and structure of the associated protein has been lacking, limiting the advances in identification of novel targets for new insecticides, especially mosquitoes. However, with more emerging data and scientific evidence in this regard in recent times, structural data of the essential proteins implicated in JH-linked pathways can be targeted and exploited. With the availability of the x-ray crystallographic structural data of mosquito juvenile hormone binding protein (MJHBP) (PDB ID: 5V13) in complex with JH ligand, this current study leverages on it and aims to identify and study in silico, potential anti-mosquito compounds targeting juvenile hormone (JH) mediated pathways via the mosquito juvenile hormone binding protein (MJHBP).

Materials and methods

Protein targets and ligands

Receptors: Crystal structure of mosquito juvenile hormone binding protein (MJHBP) (PDB ID: 5V13).

Ligands: ZINC database, African Natural Product Database (ANPDB) compounds selected using a pharmacophore model and pyrethroids.

Standalone offline softwares:

Molecular Operating Environment, 2015 version

Discovery Studio 2019 Client Full Package.

Methods

Ligand selection and preparation

Three sets of ligand were used for this study: (1) Pyrethroids, with the intention of repurposing them for control of mosquito vectors at the developmental stages, (2) Natural compounds derived from African Natural Product Database (ANPDB) (http://african-compounds.org/anpdb) and (3) compounds derived from ZINC database (http://zinc.docking.org). 3D structural data of twenty-nine (29) known pyrethroids structures were downloaded from PubChem (http://pubchem.ncbi.gov) for the study in SDF format and imported into a database in MOE. As a preliminary filter to reduce the number of compounds to be screened from ANPDB and ZINC database, a simple ligand-based pharmacophore model (Fig. 1A) was developed in MOE using the compounds depicted in Fig. 1B. The compounds were chosen based on their reported JH-target protein binding or JH-agonistic or antagonistic actions:

Fig. 1A

Pharmacophore used to filter the databases ZINC and ANPDB database B) compounds used in building the pharmacophore.

Fig. 1B

Structures of compounds used in building the pharmacophore model.

I: Juvenile hormone (a co-crystallized ligand of the target protein). It is a well known and characterized ligand of the target receptor; II and III: These are the structural derivatives of II obtained by replacing the epoxy moiety of JH with benzene and pentadiene rings respectively -the two structures showed higher binding scores relative to II on docking against MJHBP; IV and V: these two phyto-compounds isolated and purified from Lindera erythrocarpa fruits and Solidago serotina roots, and have been shown to possess strong JH-antagostic action against other JH target proteins [26]. Due to lack of information on previously identified or known inhibitors of MJHBP, the pharmacophore model was simply constructed by aligning these test compounds and selecting the common features (between two or more of the training sets) and also some unique relevant pharmacophores. The more restrictive pharmacophore model was selected. The study was limited by the non-availability of a good number of experimentally determined inhibitors and ligand thus enrichment of pharmacophore model was not possible. The study relied mostly on the co-crystallized ligand for validation. The adopted pharmacophore model (Figure A) consisted of five pharmacophoric features: Four (4) hydrophobic (Hyd) and one (1) H-bond acceptor, and was used to screen the ANPDB, and ZINC databases (Fig. 1D). The 3D - structural data of the hits were downloaded and imported into databases in Molecular Operating Environment (MOE). Compounds with very high molecular weight (>520), duplicates and huge deviations from the pharmacophore model were eliminated. Control ligands (ligand co-crystalized with receptors) were isolated from the 3D structure of the ligand-protein complex 5V13 (from RSCB PDB database) in MOE using the Seq tool. The control ligands were imported into the query ligand databases prior to preparation and docking. The database of ligand/control structures was prepared for docking as follows in MOE: protonation at a temperature of 300 K and pH 7.0 and energy minimization, using default parameters - Amber10-EHT force field was used with no periodicity, the constraints were maintained at the rigid water molecule level and partial charges were also applied [27,28].

Fig. 1D

Study design.

Protein target preparation

X-ray crystallography structural data of Mosquito Juvenile Hormone-binding protein (5V13) in pdb format was downloaded with its co-crystalized ligand (Juvenile hormone) from the RCSB PDB (Research Collaboratory for Structural Bioinformatics PDB) database (http://www.rcsb.org/pdb). The preparation and minimization was performed using tools and protocols in MOE. The preparatory process included removal of water molecules and other co-crystalized molecules, protonation, partial charges and energy minimization were implemented as described above in ligand preparation. The fully prepared and optimized 3D structure was saved in moe format for docking [[27], [28]].

Binding/docking site prediction

The co-crystalized ligand bound to the target protein defined the binding site for molecular docking. The ligand (JH) binding site option was selected as the docking site during docking simulation in MOE.

Validation of MOE docking program

A dummy docking was carried out using the x-ray crystallography structure data of the target protein (5V13) co-crystallized with a control ligand, JHIII. The isolated co-crystallized ligand was re-docked onto the binding site of the target protein. This was repeated several times with a different scoring function - ASE, Affinity dG, Alpha HB, Electron Density, GBVI/WSA dG and London dG/- each time. The docked binding pose for each scoring function was compared to the experimentally determined pose in the complex (x-ray crystallography structure). A RMSD value of less than or equal to 2.0 Å (relative to the native binding pose of the control ligand) is considered a good solution and validates a docking programme [29]. Affinity dG/GBVI/WSA dG, London dG/GBVI/WSA dG and ASE dG produced good scores and binding poses that aligned more with the experimentally determined pose of the JH-MJHBP complex (RMSD = 0.8624–1.1180 Å) (Supplementary Table 1). However, the best performance was recorded with Affinity dG scoring function (0.8624 Å) (Fig. 1C). This is also supported by evidences from Kalinowsky et al. [30]. Furthermore, the scoring function was also tested against a group of related eicosonoides – 20-hydroxyecdysone, farnesoic acid, farnesol, dinoprostone, methyl farnesoate, lipoxin A4, methroprene, JHI, JHII, and JHIII (control) - which have been shown to demonstrate relatively varied binding affinity to MJHBP in vitro (with respect to the control, JHIII) in a study by Kim et al. [18]. The scoring function correctly assigned energy score to eight (80%) out of ten of the compounds - lower negative energy score (relative to control) was assigned to five of the compounds that showed relatively little or no binding to MJHBP and about equal scores to the three that showed high binding affinity to MJHBP in Kim et al. [18]. However, two (20%) of the non-binding compounds was erroneously assigned higher negative scores by the scoring function. The Affinity dG/GBVI/WSA dG scoring function which has shown considerable good performance was used for this study. The rigid refinement approach was also adopted because the induced fit method required more computational time and did not significantly improve the predictions of the scoring functions (Supplementary Table 1).

Fig. 1C

Result of docking software validation: experimentally determined (by x-ray crystallography) pose in grey superimposed with the docked pose in green (0.8624 Å). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Docking simulation

Docking simulations were performed on MOE. The ligand was selected and docked using the Triangular matcher/rigid receptor method and scored using Affinity dG/GBVI/WSA dG options, on an Intel core i7 CPU @ 2.00 GHz, 2.60 GHz. The triangular matcher method (default in MOE) is adjudged the best placement method for standard and well-defined binding sites in MOE [31]. It generates poses by superimposing triplets of ligand atoms and triplet of receptor sites (alpha centres that represent locations of tight packing) [31]. The poses generated by the placement method was scored by the selected scoring function, Affinity dG and subsequently re-scored by GBVI/WSA dG. The Affinity dG is an empirical scoring function that calculates enthalpic contribution to binding energy using a linear function based on the following interaction factors: hydrogen bond donor−acceptor, pair, ionic interactions, metal ligation, hydrophobic interaction, and interactions between hydrophobic and polar atoms and between any two atoms [32]. The GBVI/WSA dG SF on the other hand is a force field-based scoring function which estimates the binding affinity of the ligand based on coulombic electrostatic, solvation electrostatic, van der Waals, and surface area terms trained with MMFF94x and AMBER99 force fields and ninety-nine (99) protein−ligand complexes of the solvated interaction energy (SIE) training set [33]. The protein-ligand docking poses and scores were saved in db format and ligand interaction with protein visualized (2D and 3D) using Discovery studio and MOE-ligand interaction options. The results of the top-ranking compounds were presented in tables and figures [[27], [28]].

Molecular dynamics (MD) simulation

The structures of the four highest-ranking docked protein-compound complexes (at least one from each category) were selected for a 100ns MD-simulation on a Schrodinger's Desmond module as described by Baby et al. [34]. The solvated water-soaked system was generated using the Desmond System Builder tool and the TIP3P solvating model. An orthorhombic box was simulated with a boundary distance of at least 10 Å from the outer surface of the protein with periodic boundary conditions. The system was neutralized of complex charges with the addition of 0.15 M NaCl into the simulation panel to maintain isosmotic conditions. A pre-defined equilibration procedure was performed before the simulation. The MD simulation was performed at a pressure of 1.0 bar and a temperature of 300 K (considering the target protein and organism are of insects and not human [[18], [35]], with 1000 frames saved to the trajectory, for 100 nsec period. Subsequently, the trajectory file of the simulated system was then used for calculation of various structural parameters like the Root Mean Square Deviation (RMSD), Root Mean Square Fluctuations (RMSF), Radius of Gyration (Rg), protein-ligand contacts, Intermolecular Hydrogen Bonding (H-bonding) and Solvent-Accessible Surface Area (SASA) Molecular Surface Area (MolSA), and Polar Surface Area (PSA) [34].

Binding energy calculation

The docked conformations were energy-minimized by the Primemodule of Schrodinger, and then The Molecular Mechanics Generalized Born Surface Area (MM-GBSA) analysis was implemented to calculate the binding free energy of the complexes [36].

The Formula is given below:ΔG(bind) = ΔG(solv) + ΔE(MM) + ΔG(SA)where:

ΔGsolv is the difference in GBSA solvation energy ofthe protein-ligand complex and the sum of the solvation energies for unliganded protein and ligand. ΔEMM is a difference in the minimized energiesbetween protein-ligand complex and the sum ofthe energies of the unliganded protein and ligand. ΔGSA is a difference in surface area energies of thecomplex and the sum of the surface area energies for the unliganded protein and ligand.

Result and discussion

Dengue has continued to ravage the world with a drastically growing prevalence, and incidence exacerbated by the failure of efforts and attempts to develop potent therapeutic - drug or vaccine – against the virus. The present study is a vector-targeted approach to identify potential inhibitors of juvenile-hormones action (specifically targeting mosquito juvenile hormone binding protein) as novel anti-mosquito insecticidal agents using computational techniques.

Pharmacophore modeling

The pharmacophore model filter (Fig. 1A) yielded 1930 compounds from ZINC database, and 174 from ANPDB database Fig. 1D. These compounds alongside the 29 pyrethroids were subsequently screened by molecular docking.

MOE validation

To validate our docking program, the co-crystallized MJHBP (5V13), JHIII was re-docked onto the binding site of MJHBP using different scoring functions. Affinity dG emerged as the best scoring function (Supplementary Table 1). The ligand pose/confirmation of the highest scoring pose (with Affinity dG) is presented in comparison with the experimental pose in Fig. 1C. The two poses aligned strongly with an RMSD value of 0.8624 Å. This was less than 2.0 Å and validated the program for molecular docking.

Molecular docking simulation

The study identified potential natural compounds (Table 1), synthetic compounds (Table 2) and even some pyrethroids (Table 3) with stable binding pose and very low energy score ((−10.91)- (−12.34)Kcal/mol). This was comparable to that of the control compound, Juvenile hormone III (−9.94 kcal/mol). The very low energy scores is indicative of a favourable and stabilizing binding pose of the compounds within the binding pocket of JH in the target protein, MJHBP. Studies have demonstrated that a linear relationship exists between computationally estimated binding energy and measured in vitro inhibitory activity of compounds [37,38]. The study identified a good number of compounds with high binding affinity. Binding energy alone may not be enough to ensure the inhibitory (competitive) activity of a compound. The binding pose and interaction with important active site amino acid residue to a good degree will determine the resultant effect following the binding of a compound to target protein. In the current study, we have identified compounds with extensive interaction with active site amino acids of MJHBP. The binding sites of JH III on MJHBP from x-ray crystallography studies is enclosed within a seven alpha-helices in the N-terminal domain delineated by the following amino acid residues: Leu-30, Val-48, Tyr-33, Trp-53, Tyr-64, Tyr-129, Phe-144, Val-65, Val-68, Leu-72, Leu-74, Tyr-148, and Val-51. Try-129 is reported to make a hydrogen bond interaction with the epoxy moiety of the native ligand, JH III, while Tyr-148 serves an inter-domain (N-terminal to C-terminal) linker via a hydrogen bond between its amide nitrogen and the carbonyl oxygen of Ala-285 at the C-terminus. From the result of this study, it is evident that most of the high ranking compounds from the molecular docking screening revealed interesting interactions with the amino acid residues bounding the active site of MJHBP (Table 1, Table 2, Table 3, Fig. 2, Fig. 3). Thus, the stability of their complexes with MJHBP is attributable to them: Hydrogen bonding interactions involving Trp-53 was observed commonly in most of the compounds across all three categories. For the natural compounds from the ANPDB, the hydrogen-bond interaction with MJHBP involved prominently, Trp-53, Ser-69, Tyr-129, Gly-146, and Tyr-33. The ZINC database derived compounds interaction revealed stabilizing hydrogen-bond prominently with Tyr-148 and Trp-53. Other non-covalent interactions were also observed in the complexes of these two categories of compounds. The pyrethroids formed complexes stabilized by extensive hydrophobic, alkyl and alkyl pi interactions for most of it except Flucythrinate, Tetramethrin and Imiprothrin that had hydrogen-bond interaction with amino acid residues, Tyr −148 and Trp-50.

Table 1

Molecular docking result of hit compounds from ANPDB database, on Mosquito Juvenile Hormone-binding protein (5V13).

S/N	Compound	Binding Energy	Interacting Amino Acid Residue(s)
1	Ekeberine D4	−12.3278	a(Gly146, His145), b(Gly146, Pro26), and Alkyl interactions
2	Iridin S	−12.0491	b(Trp129, Gly146, Val51,Val65, Tyr33, Tyr148, Phe144), and Alkyl interactions
3	8-methoxygenistein 7-O-alpha-l-rhamnoside 4′-O-beta glucoside	−11.6963	a(Ser69,Trp53,Tyr133,Thr29,Tyr155)Alkyl interactions
4	(E−) Pinellic acid	−11.0402	a(Tyr64,Tyr33,Tyr129,Trp53) and Alkyl interactions
5	Beta-d-glucopyranosyl 4-O-beta-d-glucopyranosyl caffeate	−11.4829	aTyr148,Ala281,Val51,Ser69,Tyr33 cTyr133, Gly146)
6	Embelin	−10.4653	a(Ser69,Val65),d(Ala281,Ala285/Val68)
7	Maesanin	−11.5640	a(Tyr129, Ser69), and Alkyl interactions
8	Subereamine A	−10.1074	a(Tyr129, Trp53,Gly146),c(Ser69),b (Trp53, Val51)

Superscripts a:Hydrogen bond, b:Carbon-hydrogen bond, c: Van der Waals.

Table 2

Molecular docking result of hit compounds from ZINC database, on Mosquito Juvenile Hormone-binding protein (5V13).

S/N	Compounds	Binding Energy (Kcal)	Interacting Amino Acid Residue(s)
1	ZINC16919139	−11.3021	a(Tyr 148, Trp 53),b(Val 51) and Alkyl interactions
2	ZINC24905365	−11.0629	a(Trp 53),b(Val 34, Ser 69, Tyr 129) and Alkyl interactions
3	ZINC16918692	−11.0048	a(Tyr 148),b(Gly 146, Val 65, Tyr 129, Pro 55, Trp 53) and Alkyl interactions
4	ZINC12944140	−10.9664	b(Val 51, Pro 26, Tyr 64, Tyr 33) and Alkyl interactions
5	ZINC16912845	−10.8944	b(Gly 146, Trp 53) and Alkyl interactions
6	ZINC16919218	−10.8237	a(Trp 53),b(Thr 29, Leu 30, Pro 26, Val 51, His 145) and Alkyl interactions
7	ZINC21523775	−10.7678	a(Tyr 148),b(Gly 146) and Alkyl interactions
8	ZINC72163779	−10.7638	b(Tyr 64, Tyr 133) and Alkyl interactions
9	ZINC00719955	−10.6818	a(Ser 69),b(Tyr 129, Ser 69, Val 65, Tyr 33) and Alkyl interactions
10	ZINC39865175	−10.6242	a(Tyr 148),b(Ala 281) and Alkyl linkages
11	ZINC93119628	−10.5189	a(Tyr 64),b(Val 51, Trp 64, Trp 50, Tyr 33, Ser 69) and Alkyl interactions
12	ZINC24906413	−10.4771	a(Tyr 148),b(Tyr 133, Val 51, Tyr 33) and Alkyl interactions
13	ZINC12886366	−10.4335	a(Trp 53),b(Trp 50, Tyr 64) and Alkyl interactions
14	ZINC24905385	−10.3880	a(Trp 50),b(Ser 69, Val 51, Val 65) and Alkyl interactions
15	ZINC16667165	−10.3421	b(Tyr 133) and Alkyl interactions
16	ZINC39863230	−10.3323	b(Trp 53, Pro 55) and Alkyl interactions
17	ZINC16919120	−10.3262	a(Trp 53),b(Leu 30, Gly 146) and Alkyl interactions
18	ZINC91742461	−10.3174	b(Tyr 33, Trp 50) and Alkyl interactions
19	ZINC24905374	−10.3133	a(Trp 53),b(Val 34) and Alkyl interactions
20	ZINC24905401	−10.2991	a(Trp 53),b(Trp 53) and Alkyl interactions
21	ZINC16919170	−10.2042	a(Trp 53),b(Gly 146, Val 65, Pro 55) and Alkyl interactions
22	ZINC39863229	−10.1999	b(Trp 53, Pro 55) and Alkyl interactions
23	ZINC24905266	−10.1854	a(Tyr 133, Trp 53),b(Trp 53, Val 65, Val 34) and Alkyl interactions
24	ZINC24906543	−10.1850	a(Tyr 148), and Alkyl interactions

Superscripts a:Hydrogen bond, b:Carbon-hydrogen bond, c: Van der Waals.

Table 3

Molecular docking result of promising pyrethroids, on Mosquito Juvenile Hormone-binding protein (5V13).

S/N	Compounds	Binding energy (S)	Interacting Amino Acid Residue(s)
1	Silafluofen	−10.9082	b(His145),d(Ile151,Tyr155,Pro26,Val51,Ala281,Trp53,/Val68,leu30) e(Tyr64, Trp53, Tyr129)
2	Esfenvalerate	−10.7287	d(Trp278,Phe269, Leu30, Val34,Phe284, Ala 285/Val51, Ala281,Ala285,Pro53) e(Tyr129,Trp53,Tyr33)
3	Flucythrinate	−10.6923	a(Tyr148),b(Val51,Tyr64),d(Ala281,Ala285,Val51/Val68), e(Tyr33, Trp53,Tyr129)
4	Resmethrin	−10.6721	b(Trp50,Val51,Tyr64),d(Tyr148,Ala285,Trp53,Phe269,Trp278,Leu30,Val34,Tyr33,Leu37,Ala281/Pro55,Val68)
5	Ethofenprox	−10.6389	d(Leu30,Pro26,Ile151,Ala251,Val51,/Ala285,Val68)e(Val51,Ala281,Ala285,Pro53)
6	Phenothrin	−10.5583	b(Trp50),d(Leu74,Phe144,Tyr129,Val68,Tyr133,Val65,Pro55,Tyr64/Ala281,Leu30)
7	Tetramethrin	−10.4103	a(Trp50),b(Trp50)d(Leu74,Phe144,Tyr133,Tyr129,Pro55,Tyr64,Val65,Val68/Val51,Phe269,Ala281,Tyr33),e(Tyr33)
8	Imiprothrin	−10.0513	a(Trp50),b(Trp50,His145)d(Tyr33,Pro55,Val65,Tyr64,Tyr129,Leu74,Val68,Phe144)
9	Control(JH)	−9.9365	a(Trp129),d(Typ278,Trp53,Ala281,Phe144,Phe269,Val51,Val48,Tyr129,Leu72,Pro55,Tyr133)

Superscripts a: Hydrogen bond, b:Carbon-hydrogen bond, c: Van der Waals, d: Alkyl/Alkyl pi Interaction, e: others.

Fig. 2

Post-docking compound-protein interactions in 3-Dimension.

Fig. 3

Post-docking compound-protein interactions displayed in 2-Dimension.

Molecular dynamics simulation

Docking studies may have been particularly successful in determining binding pose but have been shown to fail sometimes at determining ligand binding energy [39]. This is largely due to the fact that most docking protocols treat proteins as rigid molecules and assign so many approximations [40]. Thus, some docking – derived complexes may not be definite and reliable. MD simulation, which is more amenable to molecular flexibility, is used to complement docking study. In the current study, we have selected four (complexes) of the highest ranking compounds (atleast 1 from each group) from the docking results - Ekeberin D4, Maesanin, Silafluofen and ZINC16919139 (Fig. 2, Fig. 3) - for MD simulation and further analysis. Ekeberin D4 is an acyclic triterpenoid isolated from the stem-bark of the plant Ekebergia capensis [41]. It was found to possess anti-plasmodium activity and has been synthesized successfully with 27% yield in the lab [42]. Maesanin is a naturally occurring quinine isolated from species of African plants. The compound has been reported to have insecticidal activity via the inhibition the electron transport chain in flight muscle insects [43]. Silafluofen belongs to a class of insecticide, pyrethroids which have become popular as a preferred option for domestic control of insects [7]. Silafluofen was developed particularly to reduce the aquatic life-toxicity associated with pyrethroids [44]. But in recent times, resistance to this class of insecticides in adult mosquitoes has been reported [7]. ZINC16919139, on the other hand, is a synthetic compound from the ZINC database.

Stability of the protein-ligand complexes

The Root Mean Square Deviation (RMSD) and Root Mean Square fluctuations analysis is an important quantitative measure of the stability of the protein-ligand complex. A wide persistent fluctuation in RMSD and RMSF trajectory is an indicator of conformational changes in the protein/ligand system [45]. The result of this study revealed a stable system for Ekebrin D4, Silafluofen and ZINC16919139 (Fig. 4, Fig. 5). For the Ekebrine D4 complex, a major fluctuation (1.00A-2.25 Å) was observed at about 40ns in the protein (Cα), and it equilibrated henceforth with a minor fluctuations range of 1.50–2.25 Å. The ligand, on the other hand fluctuated majorly at the initial phase of the simulation, about 5ns (0.80–3.80 Å) but stabilized within a fluctuation range of 2.80–3.9 Å. Similar trends was observed with Silafluoen and ZINC16919139 - an initial phase of wide fluctuations that equilibrates and remains stable within a minor fluctuations range of <1.5 Å till the end of the simulation. ZINC16919139 ligand showed a major significant fluctuation, 1.50–4.50 Å around 40ns and equilibrated almost in that frame to a minor fluctuation range of 2.75–3.90 (1.15 Å). Silafluofen on the other hand showed a major fluctuation around 20ns (0.70–2.70 Å and 1.20–2.70 Å for both protein (Cα) and ligand respectively) Maesanin complex trajectory revealed several major fluctuations - three and four in both protein (Cα), and the ligand, respectively. There were initial major fluctuations around 18ns (1.05–2.40 Å), 30ns (1.40–2.40 Å) and 43ns (1.40–2.80 Å), which stabilized around 70ns (1.60–2.10 Å) till the end of the simulation for the protein (Cα). For the ligands, major fluctuations with wide deviations were observed very early in the simulation around 2ns, (2.2–4.2 Å), then around 30ns (2.6–4.70 Å), stabilized around 45ns (2.8–4.2 Å) and then spiked up again around the end of the simulation by (2.80–5.30 Å). Overall, the complexes seem to be fairly stable for all four compounds (all RMSD remained way less than 3.0 Å, which is considered stable). Nonetheless, Ekeberin D4, Silafluofen and ZINC16919139 appeared to have a more stable complex with MJHBP relative to Maesanin. The result of the Root Mean Square fluctuations analysis (RMSF) is consistent with that of RMSD. The RMSF value was maintained within a range of 0.40–0.25 Å for Ekeberin D4, Silafluofen and ZINC16919139 and increased to 0.50–3.70A Maesanin. This invariably implies increased perturbation of the protein system with the binding of the ligand Maesanin. The wide fluctuations in RMSF and RMSD are indications of perturbation of the system, which may be due to conformational changes within the protein complex system and repositioning of ligands inside the binding sites [40]. This has significant implication in the binding energies of the complexes (Table 4). The MMGBSA binding energies obtained for the complexes were in the following order: Ekeberin D4

Fig. 4

Post-molecular dynamics simulation analysis of protein-ligand complex (5V13-compounds) RMSD trajectory.

Fig. 5

Post-molecular dynamics simulation analysis of protein-ligand complex (5V13-compounds) RMSF trajectory.

Table 4

Binding energies (MMGBSA) of the complexes of MJHBP and the selected compound

Compounds	MMGBSA dG Bind (Kcal/mol)	MMGBSA dG Bind Coulomb (Kcal/mol)	MMGBSA dG Bind Covalent (Kcal/mol)	MMGBSA dG Bind Hbond (Kcal/mol)	MMGBSA dG Bind Lipo (Kcal/mol)	MMGBSA dG Bind Packing (Kcal/mol)	MMGBSA dG Bind Solv GB (Kcal/mol)	MMGBSA dG Bind vdW (Kcal/mol)
Silafluofen	−104.50	−7.73	−4.50E-13	−0.0003	−48.06	−4.73	23.73	−67.71
Ekeberine D4	−122.99	−31.83	−9.09E-13	−2.2203	−46.57	–	28.29	−70.67
ZINC16919139	−74.94	−8.89	−1.36E-12	−0.7037	−32.43	−1.18	32.83	−63.96
Maesanin	−72.91	−8.32	−4.50E-13	−0.5909	−34.89	–	34.29	−63.44

Radius of gyration (Rg)

The radius of gyration (Rg) is defined as the distribution of atoms of a protein around its axis. The length representing the distance between the point when it is rotating and the point where the transfer of energy has the maximum effect gives Rg. The result of this study shows that with the exception of Maesanin, the other three compounds (Ekeberin D4, Silafluofen and ZINC16919139) had a fairly stable radius of gyration - with few minor phases of fluctuation. They equilibrated with modal values of 5.30 Å, 5.85 Å, 5.57 Å, respectively (Fig. 5). Maesanin showed fluctuations around 5ns, 10ns, 30ns, and 83–90ns with a higher modal value of 6.67 Å. This finding also corroborates our earlier findings. The relatively higher radius of gyration of observed in Maesanin suggests that it has an extended structure and binding to MJHBP is associated with a conformational change that could have significantly altered the radius of gyration relative to the binding of the other compounds. This is in line with the findings of Seeliger and de Groot [46]. They reported that binding of ligand or lead compound to a protein target causes conformational changes and a consequent change in radius of gyration. An Increase in the radius of gyration (as in Maesanin) is an indication of a less compact and stable protein system [47] (see Fig. 6).

Fig. 6

Post-molecular dynamics simulation analysis of protein and ligand properties, Radius of Gyration (rGyr), Molecular Surface Area (MolSA), Solvent Accessible Surface Area (SASA) and Polar Surface Area (PSA).

The other Ligand properties, Molecular Surface Area (MolSA), Solvent Accessible Surface Area (SASA) and Polar Surface Area (PSA) are also presented in Fig. 5. They revealed consistent trajectories (with minor fluctuations) for the compounds except Maesanin. The molecular surface area (MolSA) was calculated with probe radius which was 1.4 Å and is equivalent to a van der Waals surface area [36]. The result of this study shows that MolSA plot of Ekeberin D4, Silafluofen and ZINC16919139 were fairly stable, fluctuating within the range of 485–525 Å and equilibrating at 495 Å for ekeberin; 417-427 Å and equilibrating at 423 Å for Silafluofen; 422-436 Å and equilibrating at 432 Å. Maesanin showed significant flunctuations around 10ns, 29ns, 42–56ns, 73ns and 90-96-ns. It fluctuated around range a of 390–415 Å and equilibrated at 410 Å. The Solvent Accessible Surface Area (SASA) is the surface area of a molecule accessible by a water molecule, and the Polar Surface Area (PSA) is the solvent-accessible surface area in a molecule contributed only by oxygen and nitrogen atoms. The SASA and PSA plot for this study revealed a fairly stable trajectory with few marked fluctuations observed in Ekebrin D4 around 8ns and 45–64ns PSA plot. This stabilized at 112 Å till the end of the simulation. A minor flunctuation was also observed in the early phase of PSA (0–18ns) for Silafluofen and it equilibrated around 21 Å. For ZINC16919139 and Maesanin, the equilibration was around 114A and 120 Å, respectively. The SASA plot oscillated and equilibrated for Ekeberin D4, Silafluofen Maesanin and ZINC16919139 at around 11 Å, 6.5 Å, 14.9 Å, and 3.0 Å, respectively.

Binding energy and non-covalent interactions

Non-covalent interactions (hydrophobic contacts, hydrogen bonds, and salt bridges) between proteins and ligands play are very important in stabilizing protein-ligand complex, enzyme-substrate complexes, protein-protein interaction as well nucleo-protein and significantly impact the resulting binding energy. In the present study, the MD simulation revealed that Ekeberin D4 had the highest number of H-bond interactions with relatively lesser hydrophobic interactions and water bridges (Fig. 7, Fig. 8). It made H-bond interactions with Tyr-33, Trp-50, Val-51, Trp-53, Tyr-129, Tyr-133, His-137, His-145, Gly-146, and Tyr-148. Trp-53, Tyr 129, and Tyr 133 had the highest fraction of interactions (≥60%) with the ligand. A higher fraction of hydrohphobic and water bridge interaction with Trp 64 and His 145 were observed. These non-covalent interactions may have contributed to the very high (negative) binding energy of Ekeberin D4-MJHBP complex (Table 4). In Silafluofen, we found a complex predominantly stabilized by hydrophobic interaction involving Leu-30, Tyr-33, Trp-53, Tyr-64, Val-68, Leu74, Tyr-129, Tyr-144, Tyr-268, Tyr-269, and Ala-281 - Tyr-33, Trp-53 and Val-68 predominated with over 70% fraction of interaction. This is consistent with the molecular docking result where we identified no H-bond interaction in Silafluofen-MJHBP complex. The Interaction of Ekeberin D4 with amino acid residues, Gly-146 and His-145 (observed in the molecular docking), in MJHBP were also conserved in the MD simulation. The MD simulation revealed that ZINC16919139 and Measanin showed fewer number of interactions (combination of H-bond and hydrophobic interaction) relative to Ekeberin D4 and Silafluofen. ZINC16919139 interacted with MJHBP's Trp-50, Trp53 and Tyr64 residues via H-bond (interaction fraction <60%); Leu-30, Tyr-33, Trp53, Val-68, Trp-129, Tyr-133, Ile-151 and Phe 269 via hydrophobic interaction and Val-51 and Tyr-148 via water bridge interaction. Maesanin had H-bond interaction with Tyr-33 and Tyr- 133 only (interaction fraction >50%), hydrophobic interaction with Trp53, Val-65, Tyr-129, Phe-269, Ala-285, and water bridge interaction with Tyr-33, Trp-53,Val-65, Tyr-129.

Fig. 7

The histogram of protein-ligand (5V13-compounds) contact throughout the trajectory.

Fig. 8

The total number of contacts/number of interactions in each trajectory framework of the protein –ligand complexes.

Overall, the result shows a relationship between the nature and number of interaction and the binding energy of the complexes (Table 4). Ekeberin D4-MJHBP complex with a higher number of H-bonding (and some hydrophobic) with high fraction percentage had the highest binding affinity (lowest binding energy (−122.99 kcal/mol)). While Silafluofen with extensive number of hydrophobic interaction (with high fraction percentages) equally had higher binding affinity (lower binding energy (−104.50 kcal/mol)) relative to ZINC16919139 (−74.94 kcal/mol) and Maesanin (−72.91 kcal/mol). Extensive non-covalent interactions associated with more stable complex and lower binding energy. This result is also consistent with the RMSD, and RMSF trajectory of the complexes discussed earlier. Huge fluctuations were observed in Maesanin and ZINC16919139 complexes as a result of conformational changes and ligand repositioning, (which involves loss of initial interaction and re-establishment of new ones) and this consequently impacted their overall stability and increased the binding energy.

In conclusion, this study employed an in silico structure –based approach to identify and study potential anti-mosquito Juvenile hormone binding protein (anti-MJHBP) compounds as candidates for developing anti-mosquitoe insecticides. It identified a good number of compounds with promising characteristics, high docking energy against MJHBP and capable of forming stable complexes with the protein according to MD simulation result. Further studies could prioritize and validate these compounds in in vitro studies and in vivo trials.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.