Journal of the American Chemical Society
Article
the ligand binding affinity by accounting for the effects of ligand
solvation on the protein−ligand association. The modified scoring
function to calculate the protein−ligand binding free energy in
solution (ΔGabq) can be expressed as follows.
considerable uncertainty in estimating the biochemical
potencies of putative inhibitors. Recently, X-ray crystal
structures of IKKβ have been reported both in the resting
form and in complex with a potent inhibitor.26,27 Such
structural information about the interactions between IKKβ
and small-molecule ligands can greatly aid in the design of more
effective IKKβ inhibitors.
The present study was undertaken to identify new potent
IKKβ inhibitors by consecutively applying the virtual screening
of a large chemical library to find a lead compound and de novo
design to optimize the inhibitory activity of the lead by
chemical modifications. The newly identified IKKβ inhibitors
have potential to reveal anticancer activities because the
impairment of IKKβ activity leads to the perturbation of the
inflammatory process required for cancer development. Besides
the discovery of potent IKKβ inhibitors, therefore, we also
investigated the presence of their anticancer activities with
respect to the pancreatic cancer cells.
Even when using a 3D structure of the target protein, virtual
screening and de novo design have not always been successful
due to the imperfections in the scoring function for the
estimation of the binding free energy between the target
protein and a putative ligand.28,29 One of the reasons for this
inaccuracy lies in the neglect of ligand solvation effects, which
leads to the overestimation of the binding affinity of a ligand
with many polar groups.30 Therefore, to enhance the efficiency
of virtual screening and de novo design, we modified the existing
scoring function by adding a proper molecular solvation free
energy function. The addition of this new energy term reflects
the desolvation cost for a ligand to be bound in the ATP-
binding site of IKKβ.
⎛
⎞
Aij
B
ij
ri6j ⎠
ΔGbaq = WvdW ∑ ∑ ⎜
−
⎜
⎟
⎟
12
r
ij
⎝
i=1 j=1
⎛
⎝ ij
⎞
qiqj
Cij
r12
D
ij
⎜
⎜
⎟
⎟
+ Whbond ∑ ∑ E(t)
−
+ W ∑ ∑
elec
10
ij ⎠
ε(r )r
r
ij ij
i=1 j=1
i=1 j=1
2
2
+ W Ntor + W ∑ Si(Omax − ∑ Vje−r /2σ
)
ij
tor
sol
i
i=1
j≠i
(1)
The weighting factors for van der Waals interactions, hydrogen
bonds, electrostatic interactions, and the entropic penalty, and the
molecular solvation free energy term, were represented by WvdW
,
Whbond, Welec, Wtor, and Wsol, respectively. The interatomic distance is
rij, and Aij, Bij, Cij, and Dij, are Leonard-Jones potential parameters. The
hydrogen-bond energy term includes an additional weighting factor,
E(t), to reflect the angle-dependent directionality. In calculating the
electrostatic interaction energy term, we used the sigmoidal function
ε(rij) as the distance-dependent dielectric constant of IKKβ.34 In the
torsional term, Ntor denotes the number of rotatable bonds in the
ligand. The final term in eq 1 represents the desolvation cost of the
ligand for binding to IKKβ and includes three atomic parameters, Si,
Oimax, and Vi, which denote the atomic solvation free energy per unit
volume, maximum atomic occupancy, and fragmental volume for atom
i, respectively.35 This desolvation term was calculated using the atomic
parameters reported by Park because they proved to be successful in
estimating the hydration energies of various organic molecules.36
Using this modified scoring function, we carried out docking
simulations in the ATP-binding pocket of IKKβ to score and rank
the potential IKKβ inhibitors according to the calculated binding
affinities.
Furthermore, most popular docking and de novo design
programs concentrate on finding the conformation of a ligand
that can be bound most tightly in the binding pocket of the
receptor. Although there might be a large energy difference
between this bioactive conformation of a ligand and its
minimum-energy structure,31 the conformational destabiliza-
tion of the ligand is neglected or roughly approximated from
molecular mechanics in most scoring functions. It is therefore
apparent that the inaccuracy of the scoring function can be
attributed, at least in part, to the improper description of the
ligand conformation energy. Therefore, to further enhance the
possibility of finding the potent inhibitors via de novo design, we
calculated the ligand conformational energy term in the scoring
function with ab initio quantum chemical calculations. We thus
anticipated that the present virtual screening and de novo design
procedures would enrich the chemical library with molecules
that have a good inhibitory activity against IKKβ. Finally, we
also addressed the dynamical features required for highly potent
IKKβ inhibitors based on molecular dynamics (MD)
simulations of the IKKβ−inhibitor complexes. The results
obtained in these comprehensive computational and exper-
imental studies are expected to provide useful information for
the discovery of new potent IKKβ inhibitors.
De novo Design. The structure-guided de novo design of IKKβ
inhibitors was performed in three steps. First, we made some structural
modifications to the initial IKKβ inhibitor found from the preceding
virtual screening to obtain the best scaffold from which even more
potent inhibitors could be derivatized. The LigBuilder program37 was
used in this structural transformation to obtain phenyl-(4-phenyl-
pyrimidin-2-yl)-amine (PPA) as a promising scaffold for designing new
potent inhibitors. In the second step, a variety of PPA derivatives were
generated based on the calculated binding mode of PPA in the ATP-
binding site of IKKβ. To obtain various PPA derivatives as candidates
for a potent inhibitor, a genetic algorithm was employed to generate
the derivatives by changing the chemical moieties at the substitution
positions. In this step, we used the empirical scoring function
suggested by Wang et al.38 to score and rank the generated PPA
derivatives. Because the central pyrimidin-2-yl amine group of PPA
appeared to be bound tightly to the ATP-binding site, we limited the
substitution positions to the two terminal phenyl rings. To reduce the
computational time, only the generated derivatives that could satisfy
the bioavailability rules as a drug candidate were selected for further
analysis. Based on the filtration criteria, we obtained 3617 PPA
derivatives that were estimated to have higher inhibitory activity
against IKKβ than PPA itself.
In the final de novo design step, the generated PPA derivatives were
further screened with the new scoring function constructed by
introducing a conformational destabilization energy term of ligand into
the modified scoring function shown in eq 1. This new scoring
function has the following form.
MATERIALS AND METHODS
■
Virtual Screening of IKKβ Inhibitors with Docking Simu-
lations. We used the automated AutoDock program32,33 to calculate
the binding mode and binding free energy of each molecule in the
docking library. The procedure for constructing this docking library is
detailed in Supporting Information. To improve the accuracy of virtual
screening, the original AutoDock scoring function was modified to
include a solvation free energy term for the ligand. The introduction of
this additional energy term seemed to guarantee better prediction of
ΔGbconf = ΔGbaq + ΔEconf
(2)
Here, ΔEconf represents the energy required for the conformational
change of a ligand from its ground state to the best binding mode in
the ATP-binding site of IKKβ. This new energy term was augmented
to the scoring function to enhance the accuracy of predicting the
binding free energy of a ligand and could be given by the electronic
B
J. Am. Chem. Soc. XXXX, XXX, XXX−XXX