Structure-based virtual screening for insect ecdysone receptor ligands using MM/PBSA

Article abstract of DOI:10.1016/j.bmc.2019.02.011

The ecdysone receptor (EcR) is an insect nuclear receptor that is activated by the molting hormone, 20-hydroxyecdysone. Because synthetic EcR ligands disrupt the normal growth of insects, they are attractive candidates for new insecticides. In this study, the Molecular Mechanics/Poisson–Boltzmann Surface Area (MM/PBSA) method was used to predict the binding activity of EcR ligands. Validity analyses using 40 known EcR ligands showed that the binding activity was satisfactorily predicted when the ligand conformational free energy term was introduced. Subsequently, this MM/PBSA method was applied to structure-based hierarchical virtual screening, and 12 candidate compounds were selected from a database of 3.8 million compounds. Five of these compounds were active in a cell-based competitive binding assay. The most potent compound is a simple proline derivative with low micromolar binding activity, representing a valuable lead compound for further structural optimization.

Full text of DOI:10.1016/j.bmc.2019.02.011

Accepted Manuscript
Structure-based virtual screening for insect ecdysone receptor ligands using
MM/PBSA
Shinri Horoiwa, Taiyo Yokoi, Satoru Masumoto, Saki Minami, Chiharu
Ishizuka, Hidetoshi Kishikawa, Shunsuke Ozaki, Shigeki Kitsuda, Yoshiaki
Nakagawa, Hisashi Miyagawa
PII:
S0968-0896(18)31929-1
https://doi.org/10.1016/j.bmc.2019.02.011
BMC 14747
DOI:
Reference:
Bioorganic & Medicinal Chemistry
To appear in:
Received Date:
Revised Date:
Accepted Date:
13 November 2018
1 February 2019
4 February 2019
Please cite this article as: Horoiwa, S., Yokoi, T., Masumoto, S., Minami, S., Ishizuka, C., Kishikawa, H., Ozaki,
S., Kitsuda, S., Nakagawa, Y., Miyagawa, H., Structure-based virtual screening for insect ecdysone receptor ligands
using MM/PBSA, Bioorganic & Medicinal Chemistry (2019), doi: https://doi.org/10.1016/j.bmc.2019.02.011
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers
we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and
review of the resulting proof before it is published in its final form. Please note that during the production process
errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

1
Structure-based virtual screening for insect ecdysone receptor
ligands using MM/PBSA
Shinri Horoiwa,# Taiyo Yokoi,# Satoru Masumoto, Saki Minami, Chiharu Ishizuka, Hidetoshi
Kishikawa, Shunsuke Ozaki, Shigeki Kitsuda, Yoshiaki Nakagawa,* Hisashi Miyagawa
Division of Applied Life Sciences, Graduate School of Agriculture, Kyoto University, Kyoto
606-8502, Japan
# These authors contributed equally to this work.
* To whom correspondence should be addressed.
TEL: +81-75-753-6117; FAX: +81-75-753-6123;
E-mail: naka@kais.kyoto-u.ac.jp
Abstract
The ecdysone receptor (EcR) is an insect nuclear receptor that is activated by the molting
hormone, 20-hydroxyecdysone. Because synthetic EcR ligands disrupt the normal growth of
insects, they are attractive candidates for new insecticides. In this study, the Molecular
Mechanics/Poisson–Boltzmann Surface Area (MM/PBSA) method was used to predict the
binding activity of EcR ligands. Validation study using 40 known EcR ligands showed that the
binding activity was satisfactorily predicted when the ligand conformational free energy term
was introduced. Subsequently, this MM/PBSA method was applied to structure-based
hierarchical virtual screening, and 12 candidate compounds were selected from a database of 3.8
million compounds. Five of these compounds were active in a cell-based competitive binding
assay. The most potent compound is a simple proline derivative with low micromolar binding
activity, representing a valuable lead compound for further structural optimization.
Keywords
ecdysone receptor; free energy calculation; MM/PBSA; molecular dynamics; virtual screening

2
1. Introduction
20-Hydroxyecdysone (20E; Figure 1) is an insect steroid hormone that plays a critical role in
molting.1,2 The action of 20E is mediated by the ecdysone receptor complex, a heterodimer of the
ecdysone receptor (EcR) and ultraspiracle (USP).2,3 By binding to the EcR subunit of the
complex, 20E activates the transcription of its target genes to induce molting. Because the
molting is a developmental process specific to insects, ecdysone agonists and antagonists are
attractive candidates for insecticides that are safe for mammals. In the 1980s, diacylhydrazines
(DAHs; Figure 1) were first reported as insecticidal nonsteroidal ecdysone agonists4–6 and
received considerable attention of scientists involved in insecticide discovery. To date, five DAH
analogs have been registered as insecticides to control lepidopteran and coleopteran pests in
agriculture.7 The related chemotypes, such as amidoketones (AMKs),8,9 oxadiazolines (ODZs),10
and imidazoles (IMDs),11,12 were also designed through scaffold hopping approaches (Figure 1).
However, these compounds are less potent than the commercialized DAHs and have not been
developed as insecticides.13
Figure 1. Chemical structures of representative ecdysone receptor ligands.
In 2003, Billas et al. published two X-ray crystal structures of the ligand-binding domain
(LBD) of Heliothis virescens EcR (HvEcR), one complexed with ponasterone A (PonA; Figure
1) and the other with a DAH-type ligand.14 They also crystallized another HvEcR LBD that is
bound to an IMD analog.15 These studies paved the way for structure-based rational design of
novel EcR ligands. On the basis of the active conformation of PonA at the binding site, we

3
conducted shape-based virtual screening to identify three EcR antagonists.16 Lead optimization
efforts were made in our laboratory, but the potency remained in the micromolar range
(Kitamura et al., unpublished). Recently, Hu et al. performed combined pharmacophore/docking
screening based on the three-dimensional structure of the DAH-bound HvEcR.17 The potency of
their hit compound was, however, very weak (IC50 > 100 µM). In addition, hit rates in these
studies were not so high (12.5%16 and 11.1%17). These studies underscore the need for more
effective virtual screening strategies in the EcR–ligand system.
In structure-based virtual screening, molecular docking to rigid protein models has been
widely used. However, with the rapid development of computational power, free energy
calculations based on molecular dynamics (MD) simulations, especially MM/PBSA (Molecular
Mechanics/Poisson–Boltzmann Surface Area),18 are often incorporated to increase the screening
enrichment. Several studies have demonstrated the effectiveness of the MM/PBSA method as a
virtual screening filter.19–23
The aim of this study is to validate the MM/PBSA method in the EcR–ligand system and
apply this method to large-scale structure-based virtual screening. In the preliminary step, the
binding free energy of 40 known EcR ligands was calculated using MM/PBSA and correlated
with the binding activity measured in a Spodoptera frugiperda cell line (Sf-9). In the next step,
the MM/PBSA method was incorporated into hierarchical virtual screening. From a database
containing 3.8 million compounds, 12 candidates were selected for the binding assay in Sf-9
cells. Among them, five compounds showed significant binding activity. Stereoisomers of the
most potent compound were chemically synthesized to discuss the stereochemistry–activity
relationship.
2. Materials and methods
2.1. Molecular docking and dynamics
2.1.1. Validation set
Compounds 8, 9, and 30–40 were newly synthesized as described in the Supplementary
Materials. Compounds 11–15 were generous gifts from Mitsui Chemicals Agro, Inc. (Tokyo,
Japan). The other compounds were from our stock samples.24,25 Chemical structures and binding
activities of these compounds are shown in Table 1.

4
2.1.2. Homology modeling
The three-dimensional structure of the S. frugiperda EcR (SfEcR) LBD was constructed using
MODELLER 9.16.26 Two partial sequences (GenBank accession numbers: AAM55494 and
CAD58232) were combined to yield the primary sequence of SfEcR LBD.16 The X-ray crystal
structure of the HvEcR LBD in complex with an IMD-type ligand (PDB ID: 3IXP; 2.85 Å
resolution) was used as the main template. Since the residues 311–333 are not visible in the PDB
data of 3IXP, an ecdysteroid-bound EcR structure (PDB ID: 1R1K; 2.9 Å resolution) was used
as the partial template for this region.
2.1.3. Molecular docking
Molecular docking was performed using OEDocking 3.2.0.2 software package.27 The
ligand-binding site of the SfEcR homology model was created using the MAKE_RECEPTOR
module. A conformer ensemble of each compound was generated using OMEGA 2.5.1.4,28,29
and molecular docking was performed using the HYBRID module.30,31 Docked conformations of
compounds 1–40 are shown in Figure S1. Docking models were visualized using UCSF Chimera
1.12.32
2.1.4. System preparation
AMBER16 software package33 was used for system preparation, MD simulations, and
MM/PBSA calculations. Simulation systems were built using the TLEAP module of AMBER16.
Protein–ligand complexes obtained from docking simulations were used as the initial structures.
Atomic charges for the ligands were calculated using the AM1-BCC method34,35 implemented in
the antechamber module. The AMBER force field 14SB (ff14SB)36 was applied to the protein,
and general AMBER force field 2 (GAFF2) was applied to the ligands. Na+ and Cl− ions were
added to neutralize the charge of the systems, and the systems were immersed in the cubic box of
TIP3P water molecules37 that extends at least 10 Å from the solute surface.
2.1.5. MD simulation
MD simulations were performed using the pmemd.cuda module of AMBER16. Periodic
boundary conditions were applied for all systems. Particle mesh Ewald (PME) summation38 was
employed to calculate the electrostatic interactions with a cutoff of 12.0 Å for long-range

5
interactions. The time step of the MD simulation was set to 1.0 fs. Covalent bonds connected to
hydrogen atoms were restrained by the SHAKE algorithm.39 Before the MD simulation, the
systems were energy-minimized using steepest descent and conjugate gradient algorithms to
attain a gradient of <0.5 × 10−4 kcal mol−1 Å−2. Next, the systems were gradually heated from 100
to 300 K in constant volume (NVT) ensemble over 20 ps using a restraint of 0.5 kcal mol−1 Å−2.
The systems were then equilibrated in a constant pressure (NPT) ensemble for 180 ps.
Production runs were performed in an NPT ensemble over 150 ps. Atomic coordinates of
protein–ligand complexes were stored every 1 ps for the MM/PBSA calculations. For each
ligand, MD simulation was performed in triplicate with different initial atomic velocities; thus,
450 snapshots (150 × 3) per compound were collected in total.
2.1.6. MM/PBSA calculation
As shown in Eq. 1, protein–ligand binding free energy (ΔGbind) is expressed as a difference in the
free energies of three states:
Gbind GPL GP GL
(1)
where subscripts PL, P, and L stand for protein–ligand complex, apoprotein, and free ligand,
respectively. In MM/PBSA, the free energy of each state is estimated from Eq. 2:
GX EMM GPB GSA TS
(2)
where X is PL, P, or L. In the above equation, EMM is the molecular mechanics energy in the gas
phase, which is the sum of the Coulomb interaction, the van der Waals interaction, and bond
energies. GPB is the polar term of hydration free energy, which is calculated from the
Poisson–Boltzmann equation. GSA is the nonpolar term of hydration free energy, which is
estimated from the solvent-accessible surface area. TS is the solute entropy multiplied by the
absolute temperature. Only ligand conformational free energy was considered in this study (vide
infra). The MD trajectories of the complexes were used as the source of the trajectory of the
aporeceptor and the ligand. The MM/PBSA calculation was performed using the
MMPBSA.py.MPI module of AMBER16. The interior and exterior dielectric constants were set

6
to 1.0 and 80.0, respectively.
2.1.7. Ligand conformational free energy calculation
Ligand conformational free energy (ΔGconf_ligand) was estimated using FREEFORM module of
SZYBKI 1.9.0.3.40 FREEFORM implements a rapid estimation method for ligand
conformational entropy developed by Wlodek et al.41 While the MM/PBSA calculation was
applied to 450 snapshots/compound, the conformational free energy calculation was applied to
45 snapshots/compound to reduce the computational cost.
2.2. Virtual screening
2.2.1. Database
A conformational database of 3.8 million compounds generated with OMEGA
(namiki_unfiltered_200confs.oeb.gz) was kindly provided by OpenEye Japan (Tokyo, Japan).
Each member of the database is commercially available via Namiki Shoji Co., Ltd. (Tokyo,
Japan).
2.2.2. Virtual screening
Virtual screening was performed in three steps: (i) docking screening, (ii) single-point
MM/PBSA rescoring, and (iii) MD-based MM/PBSA rescoring coupled with ligand
conformational free energy calculation (MM/PBSA + ΔGconf_ligand). Docking screening was
performed using the FRED30,31 module of OEDocking. Each member of the database was docked
into the ligand-binding pocket of the SfEcR homology model, and the docking poses were scored
by the Chemgauss4 scoring function. The top 5,000 compounds were submitted for single-point
MM/PBSA calculations. By setting a cutoff value of −10 kcal/mol, the top 389 compounds were
selected for MM/PBSA + ΔGconf_ligand rescoring. On the basis of the structural novelty and
commercial availability, 12 compounds were purchased for the binding assay. The purity of the
purchased compounds is guaranteed as >90% by the vendors. Vendor-supplied 1H NMR spectra
of the hit compounds were provided in Figure S2–S6. Structures of these compounds were
further verified by high-resolution mass spectrometry (HRMS) using a Thermo Fisher Scientific
Exactive Plus instrument (Waltham, MA, USA) equipped with an electrospray ionization (ESI)
source. Compound 41 (Asinex Ltd., Moscow, Russia), HRMS–ESI (m/z): [M + Na]+ calcd for

7
C20H21F2N3O2Na 396.1494, found 396.1492; compound 42 (Enamine Ltd., Kyiv, Ukraine),
HRMS–ESI (m/z): [M + Na]+ calcd for C21H30N2O3Na 381.2149, found 381.2150; compound 46
(Enamine), HRMS–ESI (m/z): [M + Na]+ calcd for C22H29ClN2O2Na 411.1810, found 411.1806;
compound 48 (Enamine), HRMS–ESI (m/z): [M + Na]+ calcd for C20H29N3ONa 350.2203, found
350.2204; compound 50 (Enamine), HRMS–ESI (m/z): [M + Na]+ calcd for C24H32N2O2Na
403.2356, found 403.2353. FAF-Drugs442,43 and Aggregator Advisor44 web servers were used to
detect assay interference substructures of the hit compounds.
2.3. Binding assay
The binding assay was performed as described in our earlier publications.45,46 Briefly, 400 µL of
Sf-9 cell suspension was incubated at 25 °C for 30 min in the presence of 1 µL of a test
compound in dimethyl sulfoxide (DMSO) and 2 µL of tritiated ponasterone A ([3H]PonA;
3.5–5.2 TBq/mmol, 0.5 kBq/µL) in 70% ethanol. The suspension was then filtered through a
GF-75 glass filter (Advantec Toyo Kaisha, Ltd., Tokyo, Japan) with the aid of a vacuum pump.
The filter was washed with water (3 × 3 mL), dried, and immersed in 3 mL of Insta-Gel Plus
(PerkinElmer Inc., Waltham, MA, USA) to measure the radioactivity on an LSC-6100 counter
(Aloka Co., Ltd., Tokyo, Japan). Treatments with the carrier (DMSO) and 2.7 µM unlabeled
PonA were used to determine the total and nonspecific binding, respectively. The half maximal
inhibitory concentrations (IC50, M) were calculated by probit analysis using PriProbit 1.63,47 and
the reciprocal logarithms (pIC50) were used as the index of binding activity.
2.4. Synthesis of stereoisomers of 42
2.4.1. General considerations
All reagents were purchased from commercial sources and used as received unless otherwise
stated. tert-Butoxycarbonyl (Boc)-protected L- and D-prolines [(S)- and (R)-53] were obtained
from Watanabe Chemical Industries, Ltd. (Hiroshima, Japan). (S)- and
(R)-(1,2,3,4-Tetrahydronaphthalen-1-yl)methanamines [(S)- and (R)-54] were prepared
according to the literature protcedure.48 HATU
(1-[bis(dimethylamino)methylene]-1H-1,2,3-triazolo[4,5-b]pyridinium 3-oxide
hexafluorophosphate) was purchased from Novabiochem (Burlington, MA, USA). Anhydrous
CH2Cl2 and N,N-diisopropylethylamine (DIPEA) were obtained from Nacalai Tesque, Inc.

8
(Kyoto, Japan) and Fujifilm Wako Pure Chemical Corp. (Osaka, Japan), respectively. Flash
column chromatography was performed on a Biotage Isolera One system (Uppsala, Sweden)
equipped with SNAP Ultra cartridges (silica gel, 25 µm). Melting points (Mp) were determined
on a Yanaco melting point apparatus (Kyoto, Japan) and are uncorrected. Optical rotations were
measured on a Jasco P-1010 polarimeter (Tokyo, Japan). NMR spectra was recorded on a Bruker
Avance III 400 spectrometer (Billerica, MA, USA) at 298 K. Chemical shifts (δ) are given in
ppm and calibrated using the internal standards [tetramethylsilane (0 ppm) for 1H; DMSO-d6
(39.52 ppm) for 13C]. HRMS data were obtained using the above-mentioned instrument.
2.4.2. General procedure for the HATU coupling
To a solution of 53 (0.5 mmol) and 54 (0.55 mmol) in anhydrous CH2Cl2 (1 mL) was added
HATU (0.55 mmol), followed by addition of DIPEA (1.5 mmol). The resulting mixture was
stirred at room temperature for 3 h. The mixture was diluted with CHCl3 (10 mL) and washed
successively with water (2 × 10 mL) and brine (10 mL). The organic layer was dried over
anhydrous Na2SO4 and concentrated under reduced pressure. The residue was purified by flash
column chromatography to yield the desired product.
2.4.3. tert-Butyl
(S)-2-({[(S)-1,2,3,4-tetrahydronaphthalen-1-yl]methyl}carbamoyl)pyrrolidine-1-carboxylate
[(S,S)-42]
According to the general procedure, the title compound was obtained as a colorless solid in 93%
yield (167 mg) from (S)-53 and (S)-54. Mp: 126–127 °C. [α]D18 −97.2 (c 0.21, CHCl3). 1H NMR
(400 MHz, DMSO-d6, mixture of rotamers): δ 1.34 (5.9H, s), 1.40 (3.1H, s), 1.55–1.89 (7H, m),
1.97–2.19 (1H, m), 2.59–2.78 (2H, m), 2.81–2.95 (1H, m), 3.13–3.48 (4H, m), 4.01–4.12 (1H,
m), 7.01–7.15 (3H, m), 7.15–7.23 (1H, m), 7.91 (0.35H, br t, J = 5.1 Hz), 8.05 (0.65H, br t, J =
5.6 Hz). 13C NMR (101 MHz, DMSO-d6, mixture of rotamers): δ 18.5, 18.8, 23.1, 23.9, 24.7,
28.0, 28.1, 29.0, 30.0, 31.1, 37.3, 37.5, 44.1, 44.2, 46.5, 46.6, 59.8, 59.9, 78.4, 78.5, 125.5, 125.8,
128.7, 128.9, 136.9, 137.8, 137.9, 153.4, 153.7, 172.1, 172.5. HRMS–ESI (m/z): [M + Na]+ calcd
for C21H30N2O3Na, 381.2149; found, 381.2159.
2.4.4. tert-Butyl

9
(S)-2-({[(R)-1,2,3,4-tetrahydronaphthalen-1-yl]methyl}carbamoyl)pyrrolidine-1-carboxylate
[(S,R)-42]
According to the general procedure, the title compound was obtained as a colorless gum in 94%
yield (171 mg) from (S)-53 and (R)-54. [α]D18 −72.7 (c 0.26, CHCl3). 1H NMR (400 MHz,
DMSO-d6, mixture of rotamers): δ 1.34 (5.6H, s), 1.41 (3.4H, s), 1.54–1.89 (7H, m), 1.99–2.20
(1H, m), 2.59–2.78 (2H, m), 2.82–2.94 (1H, m), 3.01–3.17 (1H, m), 3.22–3.45 (3H, m),
4.01–4.13 (1H, m), 7.01–7.15 (3H, m), 7.15–7.24 (1H, m), 7.97 (0.38H, br t, J = 5.4 Hz), 8.06
(0.62H, br t, J = 5.6 Hz). 13C NMR (101 MHz, DMSO-d6, mixture of rotamers): δ 18.5, 18.6,
23.1, 23.9, 24.6, 24.7, 28.0, 28.1, 29.0, 30.0, 31.1, 37.1, 37.3, 44.1, 44.3, 46.4, 46.6, 59.7, 59.9,
78.4, 78.5, 125.5, 125.8, 128.7, 128.9, 136.9, 137.88, 137.94, 153.4, 153.6, 172.1, 172.4.
HRMS–ESI (m/z): [M + Na]+ calcd for C21H30N2O3Na, 381.2149; found, 381.2160.
2.4.5. tert-Butyl
(R)-2-({[(S)-1,2,3,4-tetrahydronaphthalen-1-yl]methyl}carbamoyl)pyrrolidine-1-carboxylate
[(R,S)-42]
According to the general procedure, the title compound was obtained as a colorless gum in 87%
yield (161 mg) from (R)-53 and (S)-54. [α]D18 +68.4 (c 0.28, CHCl3). The 1H NMR data were in
accordance with those of (S,R)-42. HRMS–ESI (m/z): [M + Na]+ calcd for C21H30N2O3Na,
381.2149; found, 381.2160.
2.4.6. tert-Butyl
(R)-2-({[(R)-1,2,3,4-tetrahydronaphthalen-1-yl]methyl}carbamoyl)pyrrolidine-1-carboxylate
[(R,R)-42]
According to the general procedure, the title compound was obtained as a colorless solid in 84%
yield (157 mg) from (R)-53 and (R)-54. Mp: 125–126 °C. [α]D18 +92.1 (c 0.27, CHCl3). The 1H
NMR data were in accordance with those of (S,S)-42. HRMS–ESI (m/z): [M + Na]+ calcd for
C21H30N2O3Na, 381.2149; found, 381.2160.
2.5. Principal component analysis
The interaction patterns of test ligands (1–40) and top 389 compounds in the virtual screening
were generated by per-residue decomposition of their MD-based MM/PBSA scores using

10
MMPBSA.py.MPI module of AMBER16. The resulting data matrix (429 compounds × 242
residues) was subjected to principal component analysis (PCA) using scikit-learn.49

11
Table 1. Observed binding activity and calculated binding free energy of diacylhydrazines and related compounds.
Compounda)
pIC50 (M) Calculated ΔGbind (kcal/mol)
MD-based
MM/PBSA
Docking Single-point MD-based
scoref)
MM/PBSA
MM/PBSAg)
+ ΔGconf_ligand
No.
1b)
2b)
3b)
4b)
5b)
6b)
7b)
8
X
Y
R
H
H
t-Bu
6.44
7.64
7.27
6.86
6.80
7.09
8.81
8.17
8.03
7.92
7.94
8.78
8.60
6.58
6.71
7.08
6.04
6.88
-14.93
-15.15
-13.69
-13.33
-14.05
-14.03
-15.30
-16.12
-14.35
-14.26
-13.94
-15.37
-14.91
-14.21
-16.22
-13.65
-13.15
-14.91
-9.94
-8.32 (3.75)
-7.52
H
4-CH2CH3
H
t-Bu
-14.86
-13.26
-17.71
-13.84
-15.16
-16.92
-19.48
-18.71
-20.20
-19.94
-22.40
-20.78
-15.00
-16.24
-14.11
-18.79
-18.82
-10.17 (3.60) -8.94
-10.79 (3.88) -9.97
-14.72 (4.24) -13.30
3,5-(CH3)2
3,5-(CH3)2
3,5-(CH3)2
3,5-(CH3)2
3,5-(CH3)2
3,5-(CH3)2
3,5-(CH3)2
3,5-(CH3)2
3,5-(CH3)2
t-Bu
2-CH3
t-Bu
3-OH
t-Bu
-8.71 (4.11)
-7.18
3-OCH3
t-Bu
-10.37 (4.75) -8.72
-14.21 (4.10) -12.97
-16.55 (3.96) -12.17
-15.43 (4.12) -14.14
-13.91 (4.60) -12.42
-14.30 (3.47) -12.27
-15.34 (3.97) -13.42
-13.96 (3.78) -11.82
-10.65 (4.09) -9.07
-10.49 (3.69) -9.42
-10.86 (3.33) -10.07
-13.64 (3.67) -11.76
4-CH2CH3
4-CH2CH3
4-CH2CH3
2,3-(CH3)2
2-CH3-3,4-OCH2O-
2-CH3-3,4-(CH2)3O-
2-CH3-3,4-O(CH2)2O-
2-CH3-4,5-O(CH2)3-
3-CH3-4,5-O(CH2)3-
H
t-Bu
C(CH3)2CH2CH3
C(CH3)2CH(CH3)2
9
10b)
11c)
t-Bu
t-Bu
t-Bu
t-Bu
t-Bu
t-Bu
t-Bu
t-Bu
t-Bu
12 b,c) 3,5-(CH3)2
13c)
14c)
15c)
16b)
17b)
18b)
3,5-(CH3)2
3,5-(CH3)2
3,5-(CH3)2
2-Cl
2-Cl
2-NO2
2-Cl
3-CH3
-11.1 (4.42)
-9.66

12
19b)
20b)
21b)
22b)
23b)
24d)
25b)
26d)
27b)
2-Cl
4-CH2CH3
4-CN
t-Bu
t-Bu
t-Bu
t-Bu
t-Bu
t-Bu
t-Bu
t-Bu
t-Bu
7.89
5.92
7.82
5.54
8.48
8.08
7.52
6.43
6.56
-15.04
-14.44
-14.89
-15.05
-14.73
-15.41
-14.38
-12.77
-14.35
-19.39
-19.50
-20.20
-13.34
-20.11
-22.83
-18.40
-12.85
-12.73
-11.46 (3.85) -10.24
-10.74 (4.20) -9.97
-12.49 (3.82) -11.30
2-Cl
2-Cl
4-CH(CH3)2
4-SO2CH3
2-CH3-3-OCH3
2-CH3-3,4-(CH2)3O-
2,6-F2
2-Cl
-6.46 (3.81)
-5.26
2-Cl
-12.17 (3.17) -10.70
-13.21 (3.10) -11.26
-13.65 (3.53) -13.21
2-Cl
2-Cl
2-Cl-3,5-(CH3)2
3-Cl
H
-9.98 (4.08)
-9.17
H
-11.35 (3.54) -10.09
28b)
29b)
30
5.92
5.75
6.35
-13.09
-13.64
-17.35
-18.81
-18.56
-15.90
-11.16 (3.56) -9.27
-12.49 (3.50) -8.83
-13.85 (3.84) -10.56
7.22
8.12
-17.85
-18.64
-18.75
-22.08
-13.00 (4.41) -10.23
-16.16 (3.74) -12.50
31
32

13
6.70
6.49
4.88
7.27
7.95
4.50
-14.46
-16.38
-14.38
-15.76
-15.58
-12.30
-14.45
-15.54
-16.77
-18.85
-12.75
-12.26
-11.93 (3.82) -8.86
-14.38 (4.34) -8.71
-12.08 (4.42) -8.00
-14.84 (4.23) -10.85
-15.85 (4.33) -14.47
-10.30 (3.90) -6.42
33
34
35e)
36
37
38
4.80
5.05
-12.69
-12.84
-7.02
-8.16
-12.08 (3.89) -5.10
-16.01 (3.82) -9.18
39
40
a) Newly synthesized unless otherwise noted.
b) Cited from Ref.24

14
c) Kindly gifted from Mitsui Chemicals Agro.
d) Previously synthesized.25 The activity was newly evaluated.
e) The activity was measured as the racemate, but the (R)-enantiomer was used for the ΔGbind calculations.
f) Calculated using Chemgauss4 scoring function.
g) Figures in the parentheses indicate the standard deviations.

15
3. Results and discussion
3.1. Validation of MM/PBSA in EcR–ligand system
To validate the performance of MM/PBSA in the EcR–ligand system, we calculated the binding
free energy of 40 EcR ligands with various structures, namely, DAHs (1–29), IMDs (30–32),
ODZs (33, 34), and AMKs (35–40). The other chemotypes such as ecdysteroids and
tetrahydroquinolines were not included because these chemotypes are considered to have distinct
modes of binding to the receptor.14,50–52 We examined three MM/PBSA strategies (Figure
2A–C): (A) MM/PBSA score based on a single structure obtained through the MM minimization
of the complex (single-point MM/PBSA), (B) MM/PBSA score based on the MD trajectory of
the complex (MD-based MM/PBSA), and (C) the addition of the ligand conformational free
energy term to strategy B (MD-based MM/PBSA + ΔGconf_ligand). For comparison, results of the
docking scoring are also shown (Figure 2D).
A
B
-25
-20
-15
-10
-5
-25
-20
-15
-10
-5
MD-based MM/PBSA
y = −0.87x − 6.39
R2 = 0.18
Single-point MM/PBSA
y = −2.07x − 2.21
R2 = 0.38
0
0
0
2
4
6
8
10
12
0
2
4
6
8
10
12
Obsd pIC50 (M)
Obsd pIC50 (M)
C
D
-25
-20
-15
-10
-5
-25
-20
-15
-10
-5
Chemgauss4 scoring
y = −0.56x − 10.77
R2 = 0.21
MD-based MM/PBSA + ΔGconf_ligand
y = −1.50x + 0.23
R2 = 0.56
0
0
0
2
4
6
8
10
12
0
2
4
6
8
10
12
Obsd pIC50 (M)
Obsd pIC50 (M)
Figure 2. Relationships between observed pIC50 and binding free energy calculated by (A)
single-point MM/PBSA, (B) MD-based MM/PBSA, (C) MD-based MM/PBSA + ΔGconf_ligand
,

16
and (D) Chemgauss4 scoring. The solid lines represent regression lines. The R values denote the
Pearson correlation coefficients.
As shown in Figure 2B, the MD-based MM/PBSA approach shows poor consistency with the
experimental binding affinity (R2 = 0.18); surprisingly, it was even worse than the docking
scoring (Figure 2D). However, as shown in Figure 2C, introducing the ligand conformational
free energy term (ΔGconf_ligand) significantly improved the correlation (R2 = 0.56). It is noteworthy
that the slope and the intercept of the regression line are reasonable in Figure 2C. Although the
chemical structures of the test ligands are similar, the values of the conformational free energy
varied, ranging from 0.44 to 6.9 kcal/mol. This great contribution of the ΔGconf_ligand term is
probably because some ligands are forced to adopt highly unfavorable conformations at the
binding site in terms of potential energy. For example, compounds 7 and 35 have similar basic
scaffolds and adopt almost the same conformations at the binding site (Figure 3). However, the
difference in the ΔGconf_ligand values is as large as 2.84 kcal/mol, which results mainly from the
potential energy terms (see Table S1 in Supplementary Materials). This difference well explains
the fact that compound 7 is approximately 8,500 times more potent than compound 35.
Figure 3. Superposition of docked conformations of compounds 7 (tan) and 35 (light blue).
In MM/PBSA, the conformational free energy is often calculated by normal mode analysis
(NMA) based on the MD trajectory of a free ligand. However, it is challenging to apply NMA to
large-scale virtual screening due to the difficulty in convergence and the high computational cost.
In this study, FREEFORM was used to calculate the conformational free energy of the ligands.
In FREEFORM, the conformational free energy is calculated using the high-resolution
conformer ensemble generated by OMEGA28,29 and the fast estimation algorithm for ligand

17
conformational entropy developed by Wlodek et al.41 Thus, using FREEFORM does not cause
the convergence problem, and the conformational free energy of a ligand can be calculated with
an affordable computational cost for the virtual screening of a large database.
To summarize the above validation, prediction power of the examined strategies can be
ranked as follows: MD-based MM/PBSA + ΔGconf_ligand > single-point MM/PBSA > docking
scoring > MD-based MM/PBSA. On the other hand, computational efficiency can be ordered as
follows: docking scoring > single-point MM/PBSA > MD-based MM/PBSA > MD-based
MM/PBSA + ΔGconf_ligand. Therefore, we decided to use docking scoring, single-point MM/PBSA,
and MD-based MM/PBSA + ΔGconf_ligand for hierarchical virtual screening.
3.2. Virtual screening to identify novel EcR ligands
Our virtual screening strategy is outlined in Figure 4. A database of 3.8 million compounds was
screened in a hierarchical manner. First, we performed docking screening using FRED to select
the 5,000 top-ranked compounds in terms of their Chemgauss4 scores. Next, these 5,000
compounds were rescored by single-point MM/PBSA calculation, and the 389 top-scoring
compounds were chosen by setting a cutoff score of −10 kcal/mol (this value is based on the
computational resources available in our laboratory). Finally, these 389 molecules were rescored
by MD-based MM/PBSA rescoring that includes ligand conformational free energy calculation.
Compounds with the predicted ΔGbind values of ≤−9 kcal/mol were subjected to manual selection
based on their (i) structural novelty and (ii) commercial availability. Thus, 12 compounds that
are structurally distinct from the known EcR binders were selected and purchased from Namiki
Shoji Co., Ltd. (Figure 5). Although stereoisomers were treated as separate compounds in the
virtual screening, the purchased compounds are mixtures of stereoisomers.

18
Namiki library (3.8 million compounds)
Filter 1: Docking screening by FRED (Chemgauss4)
5000 compounds
Filter 2: Rescoring by single-point MM/PBSA
389 compounds
Filter 3: Rescoring by MD-based MM/PBSA + ΔGconf_ligand
12 compounds
Binding assay in Sf-9 cells
Figure 4. Outline of hierarchical virtual screening.
41
42
43
44
45
46
47
48
49
50
51
52
Figure 5. Chemical structures of compounds 41–52 selected by virtual screening. All
compounds were purchased as racemic or diastereomeric mixtures. The purity of these
compounds is >90% as reported by the vendors.

19
Table 2 summarizes the results of the competitive binding assay using [3H]PonA as the
radioactive ligand. Although seven of the screened compounds were not active even at the
highest concentrations tested (100 or 250 µM), the other compounds (41, 42, 46, 48, 50)
inhibited the radioligand binding in a concentration-dependent manner (Figure S1). The IC50
value of the most potent compound (42) was 9.1 µM, which was approximately 5,000 times
weaker than those of the commercialized compounds, tebufenozide (7) and chromafenozide (12).
However, since compound 42 is a simple proline derivative, it will serve as a versatile starting
point for further optimization.
Table 2. Calculated binding free energy and observed binding activity of screened compoundsa)
No. ΔGbind (kcal/mol)b) pIC50 (M)c)
No. ΔGbind (kcal/mol)b) pIC50 (M)c)
-9.13
-9.49
-14.01
-9.57
-9.44
-10.73
4.28
-12.99
-9.61
<3.60
4.11
41
42
43
44
45
46
47
48
49
50
51
52
5.04
<3.60
<4.00
<4.00
4.13
-10.70
-11.01
-10.16
-10.51
<3.60
4.12
<3.60
<3.60
a) All compounds were assayed as racemic or diastereomeric mixtures.
b) Calculated for a single stereoisomer using the MD-based MM/PBSA + ΔGconf_ligand protocol.
c) Average of two independent measurements.
Figure 6A shows the predicted binding pose of the most potent hit 42. Note that a single
stereoisomer of 42 [(S,S)-42] is docked into the binding site; i.e., the other stereoisomers were
filtered out through the virtual screening cascade. In this model, compound 42 is
hydrogen-bonded to the side chains of Thr57, Tyr122, and Asn218. This hydrogen-bonding
pattern is similar to that predicted for a DAH-type ligand 12 (chromafenozide; Figure 6B). The
Boc group of 42 extends toward a hydrophobic subpocket which is occupied by the tert-butyl
group of 12. Although the structure–activity relationship (SAR) of DAHs has been investigated
extensively,13 this subpocket remains poorly explored because only a few synthetic methods are
available for modifying the tert-butyl moiety of DAHs. As for compound 42, the Boc group can

20
be readily replaced with other substituents via acidic deprotection followed by acylation or
alkylation. This feature reinforces the potential of compound 42 as a lead compound.
A
B
Y122
Y122
T57
T57
N218
N218
S91
S91
Figure 6. Docking models of (A) proline-based hit 42 and (B) chromafenozide 12. The black
dashed lines indicate hydrogen bonds to Thr57, Tyr122, and Asn218.
3.3. Stereochemistry–activity relationship of proline-based hit 42
As part of the hit validation, the active compounds (41, 42, 46, 48, and 50) were inspected for the
risk of pan-assay interference compounds (PAINS)53 and colloidal aggregators54 using
FAF-Drugs442,43 and Aggregator Advisor44 web servers, respectively. All compounds were
judged neither as PAINS nor as aggregators. However, the moderate purity (>90%) of these
library compounds left open the possibility that minor impurities may account for the observed
potency. Moreover, the pIC50 values shown in Table 2 may have been underestimated by two- to
several-fold because these compounds were assayed as racemic or diastereomeric mixtures.
To resolve these ambiguities, we synthesized the four stereoisomers of the most potent hit
42 by coupling Boc-protected prolines 53 with tetrahydronaphthalene-based amines 5448 (Figure
7). The 1H NMR analyses confirmed that the synthesized stereoisomers are free of impurities
except for trace residual solvents. In DMSO-d6, each of the four stereoisomers existed as a
mixture of carbamate rotational isomers; the duplicated 1H and 13C signals were observed
(Figures S8–S13). Overlay of the 1H NMR spectra of a diastereomeric pair [(S,S)- and (S,R)-42;
Figure S14] agreed well with the vendor-supplied spectrum (Figure S3), confirming that the

21
purchased 42 is a mixture of diastereomers.
Figure 7. Chiral pool synthesis of the four stereoisomers of proline-based hit 42. RT, room
temperature.
The synthesized stereoisomers of 42 were then evaluated in the binding assay using Sf-9
cells (Table 3). Two of the four stereoisomers [(S,S)- and (R,S)-42] were twice as potent as the
purchased 42, while the other two isomers [(S,R)- and (R,R)-42] were an order of magnitude less
potent than the purchased 42. These results clearly show that (S,S)- and (R,S)-42 are responsible
for the observed potency of the purchased 42. The concentration–response curves for these
stereoisomers are shown in Figure S15.
Table 3. Calculated binding free energy and observed binding activity of stereoisomers of 42
No.
ΔGbind (kcal/mol)a) pIC50 (M)b)
(S,S)-42 -9.49
(S,R)-42 -7.15
(R,S)-42 -4.65
(S,S)-42 -3.23
5.33
4.11
5.29
4.13
a) Calculated using the MD-based MM/PBSA + ΔGconf_ligand protocol.
b) Average of two independent measurements.

22
An important aspect of the SAR data shown in Table 3 is that the configuration at the
proline α-carbon has no impact on the binding activity. Our MM/PBSA strategy failed to
rationalize this phenomenon, although the initial docked structures for the MD simulations might
be inappropriate (Figure S16). It is also possible that these stereoisomers readily epimerize under
the assay conditions, overriding the effect of the proline chiral center.
3.4. Characterization of interaction patterns of false positives
The hit rate in this virtual screening was 41.7% (five actives of 12 compounds tested), which is
markedly improved over the previous studies.16,17 However, seven compounds (43–45, 47, 49, 51,
and 52) were false positives despite their potent predicted ΔGbind values. To characterize the
false-positive behavior of these compounds, the interaction patterns of a total of 429 compounds
(40 test ligands and top 389 compounds in the screening) were analyzed using principal
component analysis (PCA). As shown in Figure 8A, most of the false positives (magenta circles)
were characterized by positive PC1 values, whereas the test ligands (green circles) and true
positives (blue circles) showed negative PC1 values. Moreover, PC1 is dominated by interaction
energy with Ser91 (Figure 8B). These results suggest that the strong attractive interaction with
Ser91 results in false positives. In fact, all false positives except for 45 were predicted to form
hydrogen bonds to the hydroxy group of Ser91 (Figures 9 and S18), while the active compounds
(1–42, 46, 48, and 50) were not (Figures 6 and S17). This serine residue is probably unavailable
for the hydrogen bonding interaction with ligands.

23
Figure 8. Principal component analysis of the ligand–receptor interaction patterns. (A) Score
plot of the first two principal components (PC1 and PC2). Green: test ligands (40 compounds),
blue: true positives (5 compounds), magenta: false positives (7 compounds), gray: top-scored
ligands in the virtual screening (377 compounds). (B) Per-residue loadings of PC1.

24
Y122
T57
N218
S91
Figure 9. Docking model of false positive 52. The black dashed lines indicate hydrogen bonds to
Thr57, Tyr122, and Asn218. The red dashed line is a hydrogen bond to Ser91, which is
characteristic of false-positive compounds.
4. Conclusions
In this study, we developed an efficient computational protocol to predict the binding activity of
EcR ligands, in which the MM/PBSA calculation based on MD trajectories is coupled to ligand
conformational free energy calculation by FREEFORM. This MM/PBSA approach was applied
to the structure-based hierarchical virtual screening, and five novel EcR ligands (IC50 < 100 µM)
were successfully identified from the database of 3.8 million compounds. The hit rate (42.7%)
was markedly improved over the previous virtual screening campaigns,16,17 highlighting the
effectiveness of our MM/PBSA approach.
Given the observed potency and ease of synthesis, the proline-based hit 42 with two chiral
centers is the most promising candidate for further optimization. To discuss the
stereochemistry–activity relationship, we prepared the four stereoisomers of 42 using the chiral
starting materials. Interestingly, the configuration at the proline α-carbon has no influence on the
binding activity. Probing the mechanism of this phenomenon, as well as optimization of the
potency, is now underway in our laboratory.
The PCA results clearly showed that false positives can be easily distinguished by their
unique interaction pattern, i.e. the hydrogen bond to Ser91. At this point, we are unable to give a
reason why this serine residue is unavailable for the hydrogen-bonding interaction with ligands.

25
Nevertheless, this information can used as a filter to remove these computational artifacts in
future virtual screening campaigns.
Acknowledgments
We would like to thank Dr. Paul Hawkins for critical reading of this manuscript; Ms. Karin
Nishimura for conducting mass spectrometric analyses; Mitsui Chemicals Agro for the gift of
chromafenozide analogs; and OpenEye Scientific Software for providing a free academic license.
Part of this study was performed in the Radioisotope Research Center of Kyoto University. This
work was supported by JSPS KAKENHI Grant Numbers JP16K07625 and JP17J01486. One of
us (T.Y.) is supported by a JSPS Research Fellowship for Young Scientists.
References
(1)
Yamanaka, N.; Rewitz, K. F.; O’Connor, M. B. Ecdysone Control of Developmental
Transitions: Lessons from Drosophila Research. Annu. Rev. Entomol. 2013, 58, 497–516.
Hill, R. J.; Billas, I. M. L.; Bonneton, F.; Graham, L. D.; Lawrence, M. C. Ecdysone
Receptors: From the Ashburner Model to Structural Biology. Annu. Rev. Entomol. 2013,
58, 251–271.
(2)
(3)
(4)
(5)
(6)
(7)
(8)
Nakagawa, Y.; Henrich, V. C. Arthropod nuclear receptors and their role in molting. FEBS
J. 2009, 276, 6128–6157.
Wing, K. D. RH 5849, a Nonsteroidal Ecdysone Agonist: Effects on a Drosophila Cell
Line. Science 1988, 241, 467–469.
Wing, K. D.; Slawecki, R. A.; Carlson, G. R. RH 5849, a Nonsteroidal Ecdysone Agonist:
Effects on Larval Lepidoptera. Science 1988, 241, 470–472.
Hsu, A. C.; Aller, H. E. Insecticidal N’-substituted-N,N’-diacylhydrazines. US 4,985,461,
1991.
Nakagawa, Y. Nonsteroidal ecdysone agonists. In Vitamins & Hormones; Litwack, G.,
Ed.; Elsevier: Amsterdam, 2005; Vol. 73, pp 131–173.
Tice, C. M.; Hormann, R. E.; Thompson, C. S.; Friz, J. L.; Cavanaugh, C. K.; Michelotti,
E. L.; Garcia, J.; Nicolas, E.; Albericio, F. Synthesis and SAR of α-acylaminoketone
ligands for control of gene expression. Bioorg. Med. Chem. Lett. 2003, 13, 475–478.

26
(9)
Tice, C. M.; Hormann, R. E.; Thompson, C. S.; Friz, J. L.; Cavanaugh, C. K.; Saggers, J. a.
Optimization of α-acylaminoketone ecdysone agonists for control of gene expression.
Bioorg. Med. Chem. Lett. 2003, 13, 1883–1886.
(10) Hormann, R. E.; Chortyk, O.; Phat Le, D. Oxadiazoline ligands for modulating the
expression of exogenous genes via an ecdysone receptor complex. US 7304162 B2, 2007.
(11) Holmwood, G.; Tietjen, K.-G.; Schindler, M.; Erdelen, C.; Lubos-Erdelen, A.;
Wachendorff-Neumann, U.; Turberg, A.; Hansen, O. Heterocyclic amide derivatives. US
2004/0248958 A1, 2004.
(12) Holmwood, G.; Schindler, M. Protein structure based rational design of ecdysone agonists.
Bioorg. Med. Chem. 2009, 17, 4064–4070.
(13) Dinan, L.; Nakagawa, Y.; Hormann, R. E. Structure–Activity Relationships of
Ecdysteroids and Non-Steroidal Ecdysone Agonists. In Advances in Insect Physiology;
Dhadialla, T. S., Ed.; Elsevier: Oxford, 2012; Vol. 43, pp 251–298.
(14) Billas, I. M. L.; Iwema, T.; Garnier, J.; Mitschler, A.; Rochel, N.; Moras, D. Structural
adaptability in the ligand-binding pocket of the ecdysone hormone receptor. Nature 2003,
426, 91–96.
(15) Moras, D.; Billas, I. M. L.; Browning, C. RCSB PDB - 3IXP: Crystal structure of the
ecdysone receptor bound to BYI08346. https://www.rcsb.org/structure/3IXP (accessed Oct
30, 2018).
(16) Harada, T.; Nakagawa, Y.; Ogura, T.; Yamada, Y.; Ohe, T.; Miyagawa, H. Virtual
Screening for Ligands of the Insect Molting Hormone Receptor. J. Chem. Inf. Model. 2011,
51, 296–305.
(17) Hu, X.; Yin, B.; Cappelle, K.; Swevers, L.; Smagghe, G.; Yang, X.; Zhang, L.
Identification of novel agonists and antagonists of the ecdysone receptor by virtual
screening. J. Mol. Graph. Model. 2018, 81, 77–85.
(18) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.;
Duan, Y.; Wang, W.; et al. Calculating Structures and Free Energies of Complex
Molecules: Combining Molecular Mechanics and Continuum Models. Acc. Chem. Res.
2000, 33, 889–897.

27
(19) Wang, J.; Kang, X.; Kuntz, I. D.; Kollman, P. A. Hierarchical Database Screenings for
HIV-1 Reverse Transcriptase Using a Pharmacophore Model, Rigid Docking, Solvation
Docking, and MM−PB/SA. J. Med. Chem. 2005, 48, 2432–2444.
(20) Kuhn, B.; Gerber, P.; Schulz-Gasch, T.; Stahl, M. Validation and Use of the MM-PBSA
Approach for Drug Discovery. J. Med. Chem. 2005, 48, 4040–4048.
(21) Okimoto, N.; Futatsugi, N.; Fuji, H.; Suenaga, A.; Morimoto, G.; Yanai, R.; Ohno, Y.;
Narumi, T.; Taiji, M. High-Performance Drug Discovery: Computational Screening by
Combining Docking and Molecular Dynamics Simulations. PLoS Comput. Biol. 2009, 5,
e1000528.
(22) Kumar, A.; Ito, A.; Hirohama, M.; Yoshida, M.; Zhang, K. Y. J. Identification of
Sumoylation Activating Enzyme 1 Inhibitors by Structure-Based Virtual Screening. J.
Chem. Inf. Model. 2013, 53, 809–820.
(23) Cao, R.; Huang, N.; Wang, Y. Evaluation and Application of MD-PB/SA in
Structure-Based Hierarchical Virtual Screening. J. Chem. Inf. Model. 2014, 54,
1987–1996.
(24) Nakagawa, Y.; Takahashi, K.; Kishikawa, H.; Ogura, T.; Minakuchi, C.; Miyagawa, H.
Classical and three-dimensional QSAR for the inhibition of [3H]ponasterone A binding by
diacylhydrazine-type ecdysone agonists to insect Sf-9 cells. Bioorg. Med. Chem. 2005, 13,
1333–1340.
(25) Soin, T.; De Geyter, E.; Mosallanejad, H.; Iga, M.; Martín, D.; Ozaki, S.; Kitsuda, S.;
Harada, T.; Miyagawa, H.; Stefanou, D.; et al. Assessment of species specificity of
moulting accelerating compounds in Lepidoptera: comparison of activity between Bombyx
mori and Spodoptera littoralis by in vitro reporter and in vivo toxicity assays. Pest Manag.
Sci. 2010, 66, 526–535.
(26) Šali, A.; Blundell, T. L. Comparative Protein Modelling by Satisfaction of Spatial
Restraints. J. Mol. Biol. 1993, 234, 779–815.
(27) OEDocking 3.2.0.2. OpenEye Scientific Software, Inc.: Santa Fe, NM, USA.
http://www.eyesopen.com (accessed Oct 30, 2018).
(28) OMEGA 2.5.1.4. OpenEye Scientific Software, Inc.: Santa Fe, NM, USA.
http://www.eyesopen.com (accessed Oct 30, 2018).

28
(29) Hawkins, P. C. D.; Skillman, A. G.; Warren, G. L.; Ellingson, B. A.; Stahl, M. T.
Conformer Generation with OMEGA: Algorithm and Validation Using High Quality
Structures from the Protein Databank and Cambridge Structural Database. J. Chem. Inf.
Model. 2010, 50, 572–584.
(30) McGann, M. FRED Pose Prediction and Virtual Screening Accuracy. J. Chem. Inf. Model.
2011, 51, 578–596.
(31) McGann, M. FRED and HYBRID docking performance on standardized datasets. J.
Comput. Aided. Mol. Des. 2012, 26, 897–906.
(32) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E.
C.; Ferrin, T. E. UCSF Chimera - A visualization system for exploratory research and
analysis. J. Comput. Chem. 2004, 25, 1605–1612.
(33) Case, D. A.; Betz, R. M.; Cerutti, D. S.; Cheatham III, T. E.; Darden, T. A.; Duke, R. E.;
Giese, T. J.; Gohlke, H.; Goetz, A. W.; Homeyer, N.; et al. AMBER 2016. University of
California: San Francisco, CA, USA.
(34) Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality
atomic charges. AM1-BCC model: I. Method. J. Comput. Chem. 2000, 21, 132–146.
(35) Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic
charges. AM1-BCC model: II. Parameterization and validation. J. Comput. Chem. 2002,
23, 1623–1641.
(36) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C.
ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from
ff99SB. J. Chem. Theory Comput. 2015, 11, 3696–3713.
(37) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L.
Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983,
79, 926–935.
(38) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth
particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577–8593.
(39) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the cartesian
equations of motion of a system with constraints: molecular dynamics of n-alkanes. J.
Comput. Phys. 1977, 23, 327–341.

29
(40) SZYBKI 1.9.0.3. OpenEye Scientific Software, Inc.: Santa Fe, NM, USA.
http://www.eyesopen.com (accessed Oct 30, 2018).
(41) Wlodek, S.; Skillman, A. G.; Nicholls, A. Ligand Entropy in Gas-Phase, Upon Solvation
and Protein Complexation. Fast Estimation with Quasi-Newton Hessian. J. Chem. Theory
Comput. 2010, 6, 2140–2152.
(42) Lagorce, D.; Sperandio, O.; Baell, J. B.; Miteva, M. A.; Villoutreix, B. O. FAF-Drugs3: a
web server for compound property calculation and chemical library design. Nucleic Acids
Res. 2015, 43, W200–W207.
(43) Lagorce, D.; Bouslama, L.; Becot, J.; Miteva, M. A.; Villoutreix, B. O. FAF-Drugs4: free
ADME-tox filtering computations for chemical biology and early stages drug discovery.
Bioinformatics 2017, 33, 3658–3660.
(44) Irwin, J. J.; Duan, D.; Torosyan, H.; Doak, A. K.; Ziebart, K. T.; Sterling, T.; Tumanian,
G.; Shoichet, B. K. An Aggregation Advisor for Ligand Discovery. J. Med. Chem. 2015,
58, 7076–7087.
(45) Nakagawa, Y.; Minakuchi, C.; Ueno, T. Inhibition of [3H]ponasterone A binding by
ecdysone agonists in the intact Sf-9 cell line. Steroids 2000, 65, 537–542.
(46) Nakagawa, Y.; Minakuchi, C.; Takahashi, K.; Ueno, T. Inhibition of [3H]ponasterone A
binding by ecdysone agonists in the intact Kc cell line. Insect Biochem. Mol. Biol. 2002,
32, 175–180.
(47) Sakuma, M. Probit analysis of preference data. Appl. Entomol. Zool. 1998, 33, 339–347.
(48) Chen, Y. K.; Bonaldi, T.; Cuomo, A.; Del Rosario, J. R.; Hosfield, D. J.; Kanouni, T.; Kao,
S.; Lai, C.; Lobo, N. A.; Matuszkiewicz, J.; et al. Design of KDM4 Inhibitors with
Antiproliferative Effects in Cancer Models. ACS Med. Chem. Lett. 2017, 8, 869–874.
(49) Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel,
M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; et al. Scikit-learn: Machine Learning in
Python. J. Mach. Learn. Res. 2011, 12, 2825–2830.
(50) Browning, C.; Martin, E.; Loch, C.; Wurtz, J.-M.; Moras, D.; Stote, R. H.; Dejaegere, A.
P.; Billas, I. M. L. Critical Role of Desolvation in the Binding of 20-Hydroxyecdysone to
the Ecdysone Receptor. J. Biol. Chem. 2007, 282, 32924–32934.