Journal of Pesticide Science
Online ISSN : 1349-0923
Print ISSN : 1348-589X
ISSN-L : 0385-1559
Regular Articles
Virtual screening identifies a novel piperazine-based insect juvenile hormone agonist
Taiyo Yokoi Taku NabeShinri HoroiwaKen’ichiro HayashiSayoko Ito-HarashimaTakashi YagiYoshiaki NakagawaHisashi Miyagawa
著者情報
ジャーナル オープンアクセス HTML
電子付録

2021 年 46 巻 1 号 p. 68-74

詳細
Abstract

Juvenile hormone (JH) agonists constitute a subclass of insect growth regulators and play important roles in insect pest management. In this work, a multi-step virtual screening program was executed to find novel JH agonists. A database of 5 million purchasable compounds was sequentially processed with three computational filters: (i) shape and chemical similarity as compared to known JH-active compounds; (ii) molecular docking simulations against a Drosophila JH receptor, methoprene-tolerant; and (iii) free energy calculation of ligand–receptor binding using a modified MM/PBSA (molecular mechanics/Poisson–Boltzmann surface area) protocol. The 11 candidates that passed the three filters were evaluated in a luciferase reporter assay, leading to the identification of a hit compound that contains a piperazine ring system (EC50=870 nM). This compound is structurally dissimilar to known JH agonists and synthetically easy to access; therefore, it is a promising starting point for further structure optimization.

Introduction

Juvenile hormones (JHs) are acyclic sesquiterpenoids that regulate various aspects of insect physiology.1) A principal role of this hormone is preventing metamorphosis induced by the molting hormone (20-hydroxyecdysone), which allows insects to proceed to larval molting. Because JHs are not present in animals other than arthropods, synthetic JH agonists (also referred to as juvenoids) are promising candidates for insecticides with reduced safety concerns. Shortly after the structure determination of JH I2) (Fig. 1), synthetic efforts by Zoecon chemists led to the discovery of the first-in-class juvenoid insecticide methoprene, which was followed by fenoxycarb (Dr. R. Maag/Roche) and pyriproxyfen (Sumitomo Chemical).3) These juvenoids were discovered through analog synthesis or random screening; no JH agonist has hitherto been developed via a rational, computational approach.

Fig. 1. Chemical structures of representative JH agonists

After the discovery of juvenoids, substantial effort has been made to elucidate their molecular mode of action. Through a mutagenesis screen in Drosophila melanogaster, Wilson et al.4,5) identified a bHLH/PAS (basic helix-loop-helix/Per-Arnt-Sim) protein termed methoprene-tolerant (Met), whose loss confers resistance to lethal doses of methoprene. Miura et al.6) found that the D. melanogaster Met (DmMet) is capable of binding JHs with nanomolar affinities; this finding was later extended to Met proteins of other insects.7,8) Furthermore, the liganded Met proteins were shown to heterodimerize with another bHLH/PAS protein, named taiman (Tai),9) and activate the transcription of early response genes, such as early trypsin10) and Krüppel homolog 1.11) These pioneering studies have established the Met/Tai complex as the intracellular receptor of JH agonists.12,13) Although the three-dimensional structure of Met remains to be solved, X-ray crystal structures of related bHLH/PAS proteins can be used to construct structural models of the JH receptor. In fact, Charles et al.7) combined homology modeling and site-directed mutagenesis to identify amino acid residues of Tribolium castaneum Met (TcMet) that are important for JH binding.

In this work, we conducted a hierarchical virtual screening program to identify novel JH agonists. This program consists of three steps: (i) ligand-based virtual screening using known JH agonists (1 and 2, Fig. 1) as the template; (ii) docking screening against the DmMet homology model; and (iii) molecular dynamics (MD)-based MM/PBSA14) (molecular mechanics/Poisson–Boltzmann surface area) rescoring corrected by the ligand conformational free energy calculation.15) From a database of 5 million compounds, 11 compounds were selected and subjected to a luciferase reporter assay in HEK293T cells that were transfected with Drosophila Met and Tai genes.16) We successfully identified a novel piperazine-based JH agonist with submicromolar potency.

Materials and methods

1. Virtual screening

1.1. Database

A structure database of 5 million commercially available compounds was provided by Namiki Shoji (Tokyo, Japan). This database was processed with OMEGA 2.5.1.417,18) to generate a maximum of 200 conformers per compound.

1.2. Shape-based screening

Shape-based screening was performed using ROCS 3.2.2.19,20) The query was generated from two known JH agonists, 2-[2-(4-phenoxyphenoxy)ethoxy]thiazole (1) and 2-[3-(4-phenoxyphenoxy)-1-propyn-1-yl)pyridine (2), using the Implicit Mills–Dean color force field. The TanimotoCombo score was calculated for each member of the database, and the top 1,000 compounds were selected for the docking screening.

1.3. Homology modeling

The three-dimensional structure of the DmMet PAS-B domain was constructed on Isolated-FAMSD system.21,22) The primary sequence of the target protein was extracted from the full-length DmMet sequence4) (GenBank accession code: AAC14350). The X-ray crystal structure of the Mus musculus Period circadian protein homolog 1 (Per1; PDB code: 4DJ2) was used as the template for the homology modeling. A Ramachandran plot of the homology model is shown in Fig. S1. The internal cavity was detected using the CASTp 3.0 web server,23) and structural figures were generated using UCSF Chimera 1.14.24)

1.4. Docking screening

Docking screening was performed using OEDocking 3.2.0.2.25) The MakeReceptor module was used to define the ligand-binding site of the Drosophila Met homology model. The FRED module26,27) was used to dock the 1,000 candidates into the binding site, and the docking poses were ranked using the Chemgauss4 scoring function. The top 120 compounds were subjected to visual inspection based on the novelty of their chemical scaffolds, and 50 compounds were selected for MD-based MM/PBSA screening.

1.5. MD-based MM/PBSA screening

MD-based MM/PBSA screening was conducted using Amber1628) and SZYBKI 1.8.0.129) according to a modified method described in the literature.15) The ligand–receptor complexes obtained from the docking screening were used as the initial structures. The atomic partial charges of the ligands were assigned using the AM1-BCC method30,31) implemented in the antechamber module of Amber16. The simulation systems were prepared using the tleap module of Amber16. The GAFF2 and ff14SB32) force fields were used to describe the ligands and the receptor, respectively. The net charge of the systems was neutralized by adding Na+ and Cl ions, and the systems were solvated in a TIP3P water box33) that extends at least 10 Å from the solute surface.

MD simulations were performed using the pmemd.cuda module of Amber16. Periodic boundary conditions were used for all simulations. Electrostatic interactions were calculated using the particle mesh Ewald method34) with a cutoff of 12.0 Å. Covalent bonds to hydrogens were restrained by the SHAKE algorithm.35) The time step was set to 1 fs. Before the MD simulations, the systems were energy-minimized using the steepest descent and conjugate gradient algorithms. The systems were then gradually heated from 0 to 300 K over 50 ps in a canonical (NVT) ensemble and equilibrated for 900 ps in an isothermal–isobaric (NPT) ensemble. Production simulations were performed for 5 ns in an NPT ensemble, during which the snapshots were sampled every 25 ps; thus, 200 snapshots were collected in each simulation.

MM/PBSA calculations were performed using the MMPBSA.py.MPI module of Amber16. The internal and external dielectric constants were set to 1.0 and 80.0, respectively. The ionic strength was set to 0.2 mM. The MM/PBSA scores were corrected by the ligand conformational free energy terms, which were calculated using the FREEFORM module36) of SZYBKI. Although the MM/PBSA calculation was performed against 200 snapshots/compound, the ligand conformational free energy calculation was applied to 20 snapshots/compound to reduce the computational cost.

2. Reporter gene assay

The reporter gene assay was performed as described in our earlier report.16) HEK293T cells were maintained in Dulbecco’s modified Eagle medium (DMEM) supplemented with 10% fetal bovine serum and 1% penicillin–streptomycin at 37°C under a 5% CO2 atmosphere. The cell suspension (400 µL/well) was seeded in a 24-well plate at a density of 0.5–1.0×105 cells/mL and incubated for 2 days. After removal of the growth medium, each well received 200 µL of the transfection mixture that contained pCS2+DmMet (100 ng), pCS2+DmTai (200 ng), pGL3-BmkJHRE×2-Fluc (180 ng), pRL-SV40 (18 ng), and the TransFast Transfection Reagent (1.5 µL; Promega, Madison, WI, USA). The cells were incubated for 1 hr, diluted with 400 µL of DMEM, and incubated for an additional 24 hr. The cells were then treated with 6 µL of a test compound solution in dimethyl sulfoxide. After incubation for 24 hr, the reporter activity was measured using the Dual-Luciferase Reporter Assay System (Promega) in accordance with the manufacturer’s protocol. The firefly luciferase activity was normalized against the Renilla luciferase activity to yield the relative luminescence unit (RLU). The half-maximal effective concentrations (EC50, M) were determined via probit regression analyses using PriProbit 1.63.37)

3. Chemical characterization of compound 13

Compound 13 was from Enamine (Kyiv, Ukraine), and its purity is guaranteed to be >90%. The structure was confirmed via 1H NMR and high-resolution mass spectrometry (HRMS) using Bruker Avance III 400 (Billerica, MA, USA) and Thermo Fisher Scientific Exactive Plus (Waltham, MA, USA) spectrometers, respectively. 1H NMR (400 MHz, 298 K, CDCl3; Fig. S2): δ 7.96–7.90 (m, 2H), 6.86–6.79 (m, 2H), 5.05–3.98 (br m, 2H), 3.87 (s, 3H), 3.82–3.10 (br m, 4H), 3.03 (br m, 1H), 1.75 (m, 1H), 1.50–1.20 (br m, 3H), 1.09–0.97 (m, 2H), 0.86–0.75 (m, 2H). HRMS–ESI (m/z): calcd for C17H23N2O3 [M+H]+, 303.1703; found, 303.1705.

Results and discussion

To identify novel JH agonists, we performed a multi-step virtual screening cascade against a database of 5 million purchasable compounds (Fig. 2). The first screening was performed using the OpenEye tool ROCS, which calculates shape and chemical (e.g., hydrogen-bond donor/acceptor, ring, hydrophobe, etc.) similarity of compounds to a query generated from known active compounds. We used compounds 1 and 2 for the query generation because these juvenoids outperformed natural JHs in our reporter system.16) The initial query generated automatically by ROCS contained six chemical features shared by these two juvenoids: i.e., three ring features (corresponding to two benzene rings and one thiazole/pyridine ring) and three hydrogen-bond acceptor features (corresponding to two ether oxygens and one thiazole/pyridine nitrogen). We removed all of the ring features and one hydrogen-bond acceptor feature to increase the structural diversity of hit compounds (Fig. 3A). All compounds in the database were ranked according to their TanimotoCombo scores, and the top 1,000 compounds were selected for the second screening.

Fig. 2. Overview of hierarchical virtual screening
Fig. 3. Structural models used for virtual screening. (A) The ROCS query generated from compounds 1 and 2. The query shape is shown as a gray transparent surface. The hydrogen-bond acceptor features are indicated by magenta mesh spheres. (B) Homology model of the DmMet PAS-B domain. The internal cavity is colored according to the Kyte–Doolittle hydrophobicity scale41) (red, hydrophobic; white, neutral; blue, hydrophilic). Cavity-lining residues are shown as thin stick models.

The second screening was done using the FRED program as the docking engine. We selected DmMet as the target for this screening because our reporter system was built on the DmMet/DmTai complex.16) Before the docking screening, we first constructed a homology model of the DmMet PAS-B domain using the X-ray structure of the mouse Per1 protein as the template. Although the sequence identity between the target and template proteins was only 20%, this was the best-matched pair at the time of modeling. The generated model contained an internal cavity surrounded mainly by hydrophobic amino acid residues (Fig. 3B), which is suitable for binding lipophilic ligands such as JHs. The 1,000 compounds that passed the first screening were docked into the cavity and ranked according to the Chemgauss4 scoring function. The top-scoring 120 compounds, whose docking scores ranged from −17 to −13, were visually inspected for the novelty of their chemical scaffolds. Thus, 50 compounds were selected for a third screening. Among these, we purchased 7 compounds (39) for assessing the reliability of the Chemgauss4 scoring.

The third screening was conducted using our MM/PBSA protocol,15) where the MM/PBSA calculation of the MD trajectory is coupled to the ligand conformer free energy estimation with FREEFORM. We performed a 5 ns MD simulation for each ligand–receptor complex selected in the second screening and extracted 200 snapshots for each complex. We then subjected these snapshots to our modified MM/PBSA calculations and estimated the free energy of ligand–receptor binding (ΔGbind). Compounds 39, which were selected based on their docking scores, showed a wide range of ΔGbind values (Table 1), and some of them were predicted to be non-binders (i.e., ΔGbind>0). We purchased four additional compounds (1013) based on their ΔGbind scores. Although stereoisomers were treated as separate compounds in the virtual screening, the purchased compounds, except for 4, were mixtures of stereoisomers.

Table 1. Docking scores, calculated binding free energy, and observed transcription-inducing activity of compounds selected in virtual screening.
No.Structurea)Chemgauss4 scorePredicted ΔGbind (kcal/mol)pEC50 (M)b)
3−16.71.69<4.00
4−16.6−1.46<4.00
5−15.8−2.24<4.00
6−15.86.48<4.00
7−15.5−0.77<4.00
8−14.93.61<4.00
9−14.13.92<4.00
10−13.9−3.79<4.00
11−13.8−5.31<4.00
12−13.3−2.57<4.00
13−13.3−2.936.06±0.11 (4)
JH I−11.8−2.578.22c)

a) All compounds except for 4 and JH I were assayed as racemic or diastereomeric mixtures. b) Mean±S.D. Values in parentheses are the number of replicates. c) Cited from Ref.16)

We submitted the 11 purchased candidates to the reporter assay in HEK293T cells16) (Table 1). The compounds were first assessed in single-dose testing at 100 µM, and only compound 13 increased the reporter activity. As shown in Fig. 4A, this compound induced clear dose-dependent activation of the luciferase gene. The EC50 value of 13 was 0.87 µM, which is 140 times less potent than JH I (Table 1). However, compound 13 (CLogP=1.28) is several orders of magnitude less hydrophobic than commercialized juvenoids (methoprene, log P=5.00; fenoxycarb, log P=4.07; pyriproxyfen, log P=5.37).38) In conjunction with the lipophilic nature of the binding site (Fig. 3B), this compound leaves much room for improving potency. In addition, the simple chemical scaffold of 13 makes optimization of its structure easy.

Fig. 4. Discovery of compound 13 as a novel JH agonist. (A) Dose-dependent activation of the luciferase reporter gene by compound 13. Each data point represents the mean±S.D. (n=3). Compound 13 was tested as a racemate. (B) Representative binding mode of 13 extracted from the MD trajectory. The hit compound, the (S)-isomer of 13, is depicted as a magenta stick model. Residues located within 4 Å from 13 are shown as thin stick models, and those hydrogen-bonded to 13 (Y404 and Y468) are highlighted in dark green. Black dotted lines indicate hydrogen-bonding interactions.

Figure 4B shows a representative docking pose of compound 13 extracted from the MD trajectory. Note that the (S)-isomer of 13 is docked into the binding site, which means that the (R)-isomer was filtered out during the virtual screening cascade. In this model, compound 13 is accommodated in the hydrophobic cavity and is hydrogen-bonded to two Tyr residues, Y404 and Y468, via its ester carbonyl and amide carbonyl oxygens, respectively. As shown in Fig. S3, these hydrogen bonds are relatively stable during the simulation, although some fluctuations were visible at the beginning. Thus, Y404 and Y468 of DmMet are likely to play crucial roles in recognizing compound 13. A sequence alignment revealed that these Tyr residues are well conserved across insects of different taxonomic orders, including Diptera, Lepidoptera, Coleoptera, Hemiptera, and Orthoptera (Fig. S4). Mutagenesis experiments against the DmMet paralog, germ cell-expressed (DmGce), revealed that Y270 in DmGce, which corresponds to Y404 in DmMet, is important for JH recognition, and docking simulations predicted a hydrogen bond between the hydroxy group of Y270 and the epoxide of JHs.39) Thus, the ester group of 13 is likely to act as a surrogate for the epoxide moiety of JHs. To test this hypothesis and improve potency, a structure–activity relationship study of 13 is currently in progress.

In this study, the 7 compounds selected based on the Chemgauss4 scores (39) were not active in the reporter assay. By contrast, one of the four compounds selected in the modified MM/PBSA rescoring (1013) significantly induced luciferase activity. Thus, our MM/PBSA protocol worked well for screening JH agonists. This is the third example where our MM/PBSA approach showed its effectiveness: the first was the virtual screening of nonsteroidal ecdysone agonists,15) and the second was the rational design of NSBR1, a nonsteroidal brassinolide-like compound.40) A major limitation of this MM/PBSA protocol is its relatively high computational cost; it cannot be used to screen millions of compounds. Nevertheless, as demonstrated in this study, incorporating this protocol into late-stage virtual screening improves the rate of screening success. We are now attempting to further improve the prediction accuracy of this MM/PBSA approach and expand its applicability, which may be reported in due course.

Conclusion

To discover novel JH agonists, we have carried out a hierarchical virtual screening campaign that consists of (i) shape-based screening using ROCS, (ii) docking screening using FRED, and (iii) modified MM/PBSA rescoring using Amber and FREEFORM. Among the 11 candidates selected, piperazine derivative 13 was successfully identified as a novel JH agonist. The submicromolar potency, moderate hydrophobicity, and simple scaffold of 13 make this compound a promising starting point for further structure optimization.

Acknowledgments

We are grateful to Dr. K. Komatsu (In-Silico Sciences) and Dr. H. Umeyama (Professor Emeritus at Kitasato University) for their help in protein modeling. We also thank OpenEye Scientific Software for providing us with an academic license. Finally, we would like to thank K. Nishimura, M. Ueno, and B. Nishikawa (Kyoto University) for supporting our spectrometric analyses. This work was supported in part by the Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Numbers JP16K07625 (to YN), JP17J01486 (to TY), and JP19K06051 (to YN). TY was supported by the JSPS Research Fellowship for Young Scientists.

Electronic supplementary materials

The online version of this article contains supplementary materials (Supplemental Figs. S1–S4), which are available at http://www.jstage.jst.go.jp/browse/jpestics/.

References
 
© Pesticide Science Society of Japan 2021. This is an open access article distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) License (https://creativecommons.org/licenses/by-nc-nd/4.0/)

This article is licensed under a Creative Commons [Attribution-NonCommercial-NoDerivatives 4.0 International] license.
https://creativecommons.org/licenses/by-nc-nd/4.0/
feedback
Top